Introduction
Head and neck cancer is a broad epithelial malignancy that arises in the paranasal sinuses, nasal cavity, oral cavity, pharynx and larynx (1). Head and neck squamous cell carcinoma (HNSCC), one type of epithelial head and neck cancer, is a heterogeneous disease and the sixth most common form of cancer worldwide (2). The 5-year survival rate of HNSCC patients with stage III or IV is extremely poor (3) and various complications will be faced by HNSCC survivors during their life time (4). Although many studies have been devoted to the exploration of methods for HNSCC diagnosis, prevention and treatment, the mechanisms of HNSCC progression remain largely unknown.
Previous studies have found that a variety of factors such as tobacco, alcohol consumption and human papillomavirus (HPV) are involved in the pathogenesis of HNSCC. For example, tobacco smoking and alcohol consumption affect the prognosis of HNSCC patients (5). HPV is involved in a subgroup of HNSCC, and is crucial for the prognosis and survival of HNSCC patients (6). In addition, increasing evidence has demonstrated that miRNAs and certain pathways play crucial roles in the etiology of HNSCC. For example, the Akt pathway (7) and Snail-RKIP signaling pathway (8) are tumorigenic targets for HNSCC therapeutic intervention. miR-21 and miR-494 have been identified as possible biomarkers for HNSCC via enhancing cell growth (9). Moreover, miR-138 has been identified as a tumor suppressor and serves as a therapeutic target for HNSCC metastasis (10). Recent studies have revealed that the prognosis of HNSCC is largely determined by tumor distribution, stage and histological characteristics at presentation. In addition, the combined treatment of cetuximab and carboplatin is promising for the survival of HNSCC patients at stage III/IV (11). Chung et al (12) identified four subclasses of HNSCC with various clinical prognoses based on microarrays. However, due to the high heterogeneous nature of these tumors, the pathogenesis of HNSCC remains unexplained by the traditional subgroups.
Microarray analysis is an effective approach to monitor global alterations in gene expression and to identify significant subtypes of HNSCC. Kim et al (13) identified the HPV status-specific significant gene set using the gene expression profile GSE39366 of HNSCC. In the present study, we used microarray analysis to screen the differentially expressed genes (DEGs) of samples in GSE3361 and GSE39366, respectively, to classify the different subtypes of heterogeneous HNSCC based on the molecular characteristics of the DEGs. Comprehensive bioinformatics was used to enrich the pathway information and miRNA-target-pathway subpath information of genes in each subtype. This study aimed to reveal the molecular mechanisms of heterogeneous HNSCC by subpath analysis and to explore several key biomarkers for the diagnosis or treatment of HNSCC in the different subtypes. This study may provide the basis for future advanced investigation into the clinical use of subtypes for HNSCC treatment.
Materials and methods
Data preprocessing and DEG screening
The gene expression profile of GSE6631 (14), containing 44 paired (from the same patient) samples of HNSCC and normal tissues, was downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database in the National Center for Biotechnology Information (NCBI) (15) based on the platform of GPL8300 Affymetrix Human Genome U95 ver. 2 array. Patients who received previous treatment (radiotherapy or chemotherapy) for the index tumor or another head and neck primary tumor within the past 5 years were excluded. Another gene expression profile of GSE39366 (16), containing 168 HNSCC samples, was downloaded from the GEO database in NCBI based on the platform of GPL9053 Agilent-UNC-custom-4X44K. The clinical characteristics of the patients included in the present study represent a broad cross-section of patients with HNSCC that is highly representative. Moreover, there was no correlation of tumor subtype with age, gender, race, alcohol use, pack years of smoking or tumor size.
All the CEL files obtained from the two gene expression profiles were preprocessed using the Python procedure (17) and then transformed into gene symbols. The mean expression value was considered as the gene expression value of each gene. The DIDS (Detection of Imbalanced Differential Signal) algorithm (18) was used to screen the DEGs of HNSCC in GSE6631 with P<0.05.
In addtion, Limma package in Bioconductor (19) was used to select the DEGs of the HNSCC samples in GSE6631 with P<0.05. DIDS algorithm was also used to screen the DEGs of HNSCC in GSE39366 with P<0.05.
Identification of HNSCC-associated DEGs and hierarchical clustering analysis
All selected genes from the GSE39366 profile were matched with the screened DEGs using the DIDS algorithm in GSE6631 to select the common DEGs in the two microarrays. The common DEGs were defined as the DEGs that were associated with HNSCC and to distinguish the 138 HNSCC samples into different subtypes based on their molecular characteristics.
In addition, Cluster software (20) was used to cluster the 138 HNSCC samples based on the gene expression values of the selected common DEGs. In addition, the samples were subtyped based on clustering analysis. Heat maps between the gene expression values and samples were generated using TreeView (21).
DEGs in different subtypes
In order to identify the significant DEGs among the different subtypes, we distributed the DEGs in each subtype based on their mean expression values using the following steps: g stands for one gene and score stands for the gene expression value of g in DIDS. If the score was >0, then g was upregulated in the case sample, otherwise, g was downregulated (score <0). The mean expression values for g in m subtypes are shown as {x1, x2,…..xi……xm}. Max and min stand for the maximum and minimum expression values of g, respectively. The interval = max − min. The following formula 1 was used to identify whether g was a specific gene in subtype i or not:
U
=
max
−
min
(
U
>
0.02
)
{
x
i
>
max
−
0.25
U
(
if score
>
0
)
x
i
<
min
+
0.25
U
(
if score
<
0
)
Whereas U stands for the distance between the max and min of g in m subtypes, a larger u represents the significance of g among the m subtypes. One gene with this type of U was determined as the specific gene of the significant subtype. Otherwise, one gene was determined as the common gene for the m subtypes.
In addition, the U distribution curve was used to identify the specific genes and common genes (Fig. 1). The max u of the genes was 0.04, while the min U was 0. Genes with U>0.02 were determined to be significant specific genes in each subtype, while genes with U<0.02 were defined as common genes in the common gene set.
Significant pathways of the specific DEGs
The molecular signatures database (MSigDB) is a collection of annotated gene sets for use with GSEA software (22). The Kyoto Encyclopedia of Genes and Genomes (KEGG) symbols and 186 KEGG pathways were downloaded from the MSigDB database to enrich the significant pathways of the specific DEGs in each subtype using the binomial distribution (23) which were suitable for enrichment analysis of a large scale of genes. The following formula 2 is shown as:
![]()
![]()
If there were M genes in one KEGG pathway, then Mi genes were the specific genes. Whereas P stands for the percentage of Mi-specific genes in the total M genes, and Pe is the background frequency; N is the background gene sets (all genes in the gene expression profile), and N1 is the number of specific genes. P<0.05 was chosen as the threshold.
Identification of the miRNA-target-pathway subpath
The abnormal expression of genes leads to abnormal pathways in disease; however, these abnormal DEGs are regulated by various miRNAs. We analyzed the relationship between important miRNAs from the miRNA database and target genes from each subtype or common gene set to select the miRNA-target-pathway subpath (Fig. 2).
The followed algorithm 3 was used to identify the important subpaths:
![]()
![]()
![]()
![]()
Score
=
log
weight
1
*
weight
2
*
weight
3
Whereas weight1 is the weight for miRNA, G¿ is the total number of specific genes, G is the number of specific genes that are regulated by miRNA. If specific genes were not regulated by miRNAs, then weight1 = 1. Also, a higher weight1 stands for the close regulatory relationship between miRNA and specific genes.
Weight2 is the weight for target gene, P¿ is the total number of pathways that all the targets are involved in, and P is the number of pathways that one target is involved in. If the target genes were not involved in any pathway, then weight2 = 1. Moreover, a higher weight2 indicates numerous pathways that one target is involved in.
Weight3 stands for the weight for the pathway, M¿ is the total number of genes that are involved in this pathway, M is the number of specific genes in this pathway. If there was no specific gene in one pathway, then weight3 = 1. Moreover, a higher weight3 represents a significant correlation between the pathway and subtype.
In addition, weight is the weight of the miRNA-target-pathway, and the score represents the weight for one subpath. A higher score stands for the significant correlation between subpath and subtype. Finally, a subpath where one miRNA regulates several genes and one gene is involved in several pathways was recognized as a significant subpath for HNSCC.
Results
DEG screening
In total, 2,528 DEGs from samples in GSE6631 were screened using DIDS algorithm with P<0.05 (Fig. 3). Moreover, 1,095 DEGs from the HNSCC samples in GSE6631 were screened compared with the normal samples using the traditional Limma in Bioconductor with P<0.05 (Fig. 3). In addition, the common DEGs selected using the two methods are shown using a Venn plot. The results showed that 1,032 genes were the common DEGs selected using the two different methods. The 63 genes screened by the Limma method were recognized as DEGs while they were non-DEGs using the DIDS method. The other 1,496 genes were the DEGs selected by DIDS method while they were non-DEGs using Limma, indicating that 94% of the DEGs in the different groups that were screened using Limma were identified with DIDS. However, the DEGs in one group screened by Limma were not able to be identified by DIDS.
Identification of HNSCC-associated DEGs and hierarchical clustering analysis
In total, 2,443 DEGs from the 138 samples in GSE39366 were screened by matching with the 2,528 DEGs selected using the DIDS algorithm from samples in GSE6631. From the heat maps, the 138 HNSCC samples were clustered into 4 subtypes based on the molecular characteristics of the selected 2,443 DEGs (Fig. 4). In addition, there were 14 samples characterized by HPV and CROI (conference on retroviruses and opportunistic infections). Thus, 10 samples were distributed in subtype 2, while the other 4 were distributed in subtype 1 and 3. Enhance, we speculated that HPV infection might be associated with HNSCC and that it may participant in the pathogenesis of subtype 2.
DEG distribution and pathway enrichment analysis
Based on formula 1, 377 specific DEGs in subtype 1, 46 specific DEGs in subtype 2, 471 specific DEGs in subtype 3 and 451 specific DEGs in subtype 4 were identified. The other 1,539 DEGs, which were not enriched in any subtype acted as the common gene set.
In addition, the enriched significant pathways of DEGs in four subtypes and in the common gene set are shown in Table I. The specific DEGs in subtype 1 were mainly enriched in the metabolism and cancer pathways, such as fatty acid metabolism, melanoma and colorectal cancer (Table IA). Specific DEGs in subtype 2 were mainly enriched in the immune related pathways, such as endocytosis, autoimmune thyroid disease, and antigen processing and presentation (Table IB). Moreover, significant pathways of specific DEGs in subtype 3 included the metabolism pathway and various cancer-related pathways, such as focal adhesion, pathways in cancer, extracellular matrix (ECM) receptor interaction, and ERBB signaling pathway (Table IC). In addiiton, the specific DEGs in subtype 4 mainly function in endocytosis, allograft rejection, cytochrome P450 and leukocyte transendothelial migration pathway (Table ID). In addition, the enriched pathways of DEGs in the common gene set were DNA replication, cell cycle and oocyte meiosis (Table IE).
Identification of significant miRNA-target-pathway subpaths
To identify the important miRNAs that are associated with the selected specific DEGs in the 4 subtypes and with DEGs in the common gene set, the miRNA-target-pathway subpaths with the top 10 scores in the 4 subtypes of HNSCC and in the common gene set were analyzed (Table II). The results showed that miRNA-target-pathway subpaths with the highest score in each subtype or in the common gene set were coincidence with the enriched significant pathways based on the binomial distribution. The top 1 miRNA-target-pathway subpath in subtype 1 was primary bile acid biosynthesis (Table IIA), the top 1 subpath of subtype 1 was sulfur metabolism (Table IIB), and the top 1 subpath of DEGs in subtype 3 was primary bile acid biosynthesis (Table IIC), while the top 1 subpath in subtype 4 was the ECM receptor interaction pathway (Table IID). In addition, the top 1 subpath in the common gene set was DNA replication and base excision repair (Table IIE). Each miRNA in one subpath regulates several target genes, indicating that the miRNA might be a potential biomarker for HNSCC. Furthermore, target genes participate in a variety of pathways, implying that they may be therapeutic targets.
Discussion
HNSCC is a heterogeneous disease and is the sixth most common form of cancer worldwide (2). Subtype analysis of significant miRNAs and genes associated with HNSCC will be of great significance for individualized clinical treatment and prevention of HNSCC. In the present study, we classified four subtypes of HNSCC and analyzed the subpaths of each subtype of HNSCC based on the regulatory relationship between miRNAs and biological pathways.
Our results revealed that 10 samples of patients infected with HPV were distributed in subtype 2 of HNSCC, indicating that HPV may be associated with HNSCC. HPV has been identified as an etiologic agent for oropharyngeal carcinoma, which is a subset of HNSCC (24). In addition, the majority of HNSCC patients have the potential to be infected with HPV, indicating a significant correlation between HPV and HNSCC (25). Moreover, HPV-positive status affects the therapeutic response and survival of HNSCC patients (26). Based on our data, we suggest that HPV infection may be associated with HNSCC and may increase the risk of HNSCC, which may provide a strategy for HNSCC prevention.
Our study showed that the miRLet-7A-CYP46A1-primary bile acid biosynthesis pathway was the top 1 subpath in subtype 1 and subtype 3 of HNSCC (Table IIA and C). CYP46A1 (cholesterol-24S-hydroxylase) is a CYP450 superfamily enzyme that can convert cholesterol to 24S-hydroxycholesterol (27). Bile acids are derived from cholesterol and synthesis of bile acid is the predominant metabolic pathway for catabolism (28), and certain bile acid metabolism selectively increases cancer risk (29). CYP7A1 (the homologue of CYP46A1) was found to be regulated by the bile acid-activated JNK pathway in primary rat hepatocytes (30). Thus, we speculate that CYP46A1 may function in HNSCC formation via the primary bile acid biosynthesis pathway. On the other hand, Let-7A is a novel biomarker for HNSCC (31). Downregulation of Let-7A has been found in many cancers (32), and it has been reported as a tumor suppressor in head and neck cancer via eliminating putative tumor-initiating cells (33). Based on our results, we speculate that miRLet-7A is a tumor suppressor for HNSCC, and regulates the target CYP46A1 via the primary bile acid biosynthesis pathway in subtypes 1 and 3 of HNSCC.
In the present study, our results revealed that the miR-1/miR-153/miR-206-BPNT1-sulfur metabolism pathway was the top 1 subpath in subtype 2 of HNSCC (Table IIB). Downregulation of miR-206 was previously found to contribute to laryngeal cancer proliferation and invasion via regulation of VEGF expression (34), while VEGF was reported to be a predictor for HNSCC (35), implying that miR-206 may be an inhibitor for HNSCC. A previous study found that miR-1 inhibited cell proliferation of HNSCC by targeting TAGLN2 (36), and downregulation of miR-153 promoted tumor metastasis in human epithelial cancer by targeting ZEB2 (37). On the other hand, nucleotide biosynthesis is involved with tumor cell growth during cell proliferation (38). Moreover, nucleotide biosynthesis and protein are coupled by PRPS2 in diving tumorigenesis (39), suggesting that nucleotide biosynthesis may be involved in HNSCC progression. Bisphosphate 3′-nucleotidase 1 (BPNT1) is a member of the magnesium-dependent phosphonoesterase family (40). The roles of BPNT1 in HNSCC have not been fully elucidated. However, BPNT1 plays a role in nucleotide metabolism (41). Hence, we infer that miR-1/miR-153/miR-206 may be promoters in subtype 2 HNSC oncogenesis and metastasis by targeting BPNT1 via the sulfur metabolism pathway.
The miR-506-COL5A1-ECM receptor interaction pathway was found to be the top 1 subpath in subtype 4 HNSCC (Table IID). Studies have revealed that miR-506 affects the EMT process by regulating various EMT-related genes in cancers, such as breast and ovarian cancer (42,43). In additon, EMT is a functional property of HNSCC metastasis (44). Therefore, miR-506 may be a contributor to HNSCC metastasis. A previous study revealed that tumor cells interact with the ECM during tumor invasion or metastasis (45). COL5A1 (collagen type VI) is a collagen family protein that is involved with ECM formation (46). Overexpression of collagen type IV is found to be positively correlated with tumor size in breast cancer (47). Based on our results, we deduce that miR-506 may play crucial roles in the promotion of HNSCC metastasis by regulating COL5A1 via the ECM receptor interaction pathway.
In addition, the miR-519A-MCM7-DNA replication pathway was found to be the top 1 subpath in the common gene set (Table IIE). Studies revealed that miR-519A was downregulated by ΔNp63α in cancers (48), and overexpres-sion of ΔNp63α was beneficial for HNSCC resistance to chemotherapy-induced cell death via inhibiting the transcription of TAp73b (49). One study found that the downregulation of ΔNp63α was a determinant of cellular response to DNA damage in HNSCC (50). Hence, miR-519A may be a tumor suppressor for HNSCC. Abnormal DNA replication promotes the instability of the genome during the early, middle or late S phase (51), while genomic instability plays crucial roles in cancer susceptibility (52), suggesting that DNA replication proteins are involved in HNSCC carcinogenesis. Minichromosome maintenance complex component 7 (MCM7) is a protein belonging to the highly conserved minichromosome maintenance family that are related to DNA replication (53). MCM2 (a member of MCM family) expression is found to be aberrant in laryngeal squamous epithelial lesions (54). In addition, MCM7 is identified as a useful biomarker for HPV-positive HNSCC (55). Our data indicate that miR-519A may act as a therapeutic biomarker for HNSCC by targeting the DNA replication associated protein MCM7.
In summary, our study revealed that there were four subtypes of HNSCC. Several miRNAs were identified as biomarkers or promoters for HNSCC in the different subtypes, followed by their relevant gene targets. miRLet-7A was found to be involved in subtype 1 and subtype 3 by targeting CYP46A1 via the primary bile acid biosynthesis pathway. miR-1/miR-153/miR-206 was found to be involved in subtype 2 with target BPNT1 via the sulfur metabolism pathway, and miR-506 in subtype 4 with COL5A1 via the ECM receptor interaction pathway. In addition, miR-519A was found in the common gene set for HNSCC with MCM7 via the DNA replication pathway. Moreover, HPV may be positively correlated with subtype 2 HNSCC and may increase the risk of HNSCC. Our study may help to explain the specific mechanisms of the HNSCC subtypes and may provide a basis for future investigation of clinical HNSCC treatments. However, further experimental studies are still needed to confirm our results.