Submit an Article
Online ISSN: 1945-4589
Impact Journals, LLC. is the publisher of Aging: www.impactjournals.com.
Impact Journals is a member of the Wellcome Trust List of Compliant Publishers.
Impact Journals is a member of the Society for Scholarly Publishing.
Learn about our FREE Post-Publication Promotion Services
Research Paper Volume 12, Issue 1 pp 767—783

Six-gene signature for predicting survival in patients with head and neck squamous cell carcinoma

Juncheng Wang1, *, , Xun Chen2, *, , Yuxi Tian3, , Gangcai Zhu4, , Yuexiang Qin5, , Xuan Chen6, , Leiming Pi7, , Ming Wei1, , Guancheng Liu8, , Zhexuan Li1, , Changhan Chen1, , Yunxia Lv9, , Gengming Cai10, ,

  • 1 Department of Otolaryngology, Head and Neck Surgery, Xiangya Hospital, Central South University, Changsha 410008, People’s Republic of China
  • 2 Department of Oral and Maxillofacial Surgery, First Affiliated Hospital of Quanzhou, Fujian Medical University, Quanzhou 362000, People’s Republic of China
  • 3 Department of Oncology, Xiangya Hospital, Central South University, Changsha, 410008 People’s Republic of China
  • 4 Department of Otolaryngology, Head and Neck Surgery, The Second Xiangya Hospital, Central South University, Changsha 410011, People’s Republic of China
  • 5 Department of Health Management, The Third Xiangya Hospital, Central South University, Changsha 410011, People’s Republic of China
  • 6 Department of Stomatology, Changzheng Hospital, Second Military Medcial University, Shanghai 200003, People’s Republic of China
  • 7 Department of Otolaryngology, Head and Neck Surgery, HeYuan People's Hospital, Jinan University, He Yuan 517000, People’s Republic of China
  • 8 Department of Otolaryngology, Head and Neck Surgery, Affiliated Hospital of Guilin University, Guilin 541000, People’s Republic of China
  • 9 Department of Thyroid Surgery, The Second Affiliated Hospital of Nanchang University, Nanchang 330006, People’s Republic of China
  • 10 Department of Otolaryngology, Head and Neck Surgery, First Affiliated Hospital of Quanzhou, Fujian Medical University, Quanzhou 362000, People’s Republic of China
* Equal contribution

Received: September 16, 2019       Accepted: December 24, 2019       Published: January 12, 2020      

https://doi.org/10.18632/aging.102655
How to Cite

Copyright: © 2020 Wang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

The prognosis of head and neck squamous cell carcinoma (HNSCC) patients remains poor. High-throughput sequencing data have laid a solid foundation for identifying genes related to cancer prognosis, but a gene marker is needed to predict clinical outcomes in HNSCC. In our study, we downloaded RNA Seq, single nucleotide polymorphism, copy number variation, and clinical follow-up data from TCGA. The samples were randomly divided into training and test. In the training set, we screened genes and used random forests for feature selection. Gene-related prognostic models were established and validated in a test set and GEO verification set. Six genes (PEX11A, NLRP2, SERPINE1, UPK, CTTN, D2HGDH) were ultimately obtained through random forest feature selection. Cox regression analysis confirmed the 6-gene signature is an independent prognostic factor in HNSCC patients. This signature effectively stratified samples in the training, test, and external verification sets (P < 0.01). The 5-year survival AUC in the training and verification sets was greater than 0.74. Thus, we have constructed a 6-gene signature as a new prognostic marker for predicting survival of HNSCC patients.

Introduction

Head and neck cancer is a highly heterogeneous malignant tumor, which can originate from various anatomical sites in the upper airway and digestive tract, including the mouth, larynx and pharynx. Most cases of head and neck cancer are squamous cell carcinoma (HNSCC), which accounts for about 4% of all new cancer diagnoses in the United States. Worldwide, there are about 600,000 new head and neck cancer patients each year [1], and the 5-year survival rate is only 40-50% [2]. The main reason for this high mortality is the high rate of diagnosis of advanced cancers, as the survival rate among advanced cancer patients is only 34.9% [3]. There is thus an urgent need for markers to help clinicians make accurate early HNSCC diagnoses, predict clinical outcomes, and provide reference for individualized medicine.

Many studies have been carried out in an effort to find predictive biomarkers to establish guidelines for the long-term prognosis of HNSCC patients. These biomarkers can be divided into two categories: 1) single molecules such as squamous cell carcinoma antigen (SCC-A), human papilloma virus (HPV), or any of the other new markers currently being studied; and 2) gene markers identified by analyzing high-throughput gene expression profiles and constructed using from several to dozens of prognostic genes. Several systems biology methods are currently being used to identify gene biomarkers related to the HNSCC prognosis and to characterize those genes [4–6]. For example, Tian et al. used weighted gene correlation network analysis and least absolute shrinkage and selection operator Cox regression to identify a 6-lncRNA signature within a gene expression profile. De et al. used gene expression meta-analysis to identify a 172-gene signature. And Zhao et al. used protein-protein interaction network analysis and Cox regression analysis to identify a 3-gene signature. All three authors tested their genetic signature using independent external data sets. However, the AUC for Zhao et al's genetic signature is not high (AUC=0.6), which means that identifying robust lncRNA signatures is still a challenge, and more investigation will be required to verify signatures. In other words, there is still an important need to identify new gene signals related to the prognosis of HNSCC through bioinformatics analysis of their biological functions.

To effectively identify a reliable gene signature associated with prognosis in HNSCC, we proposed a systematic pipeline to identify HNSCC-related gene markers. This approach enabled us to identify a 6-gene signature that can be used to effectively predict prognostic risk in HNSCC patients and provide a basis for better understanding of the molecular mechanisms underlying the prognoses of HNSCC patients.

Results

Identification of gene sets associated with total patient survival

For The Cancer Genome Atlas (TCGA) training set samples, we used single factor regression analysis to establish the relationship between overall survival (OS) and gene expression. We identified 425 single factor Cox regression logrank P values less than 0.01. Of those, 141 genes were associated a hazard ratio (HR) > 1, and 284 with a HR < 1. The 20 genes with the highest HRs are listed in Table 1.

Table 1. List of the most relevant 20 genes.

ENSG IDSYMBOLHRcoefficientz scoreP
ENSG00000127084TLL20.630-0.462-4.5625.07E-06
ENSG00000254656FAM69A1.3600.3084.2022.64E-05
ENSG00000041515HIST2H3PS21.3370.2904.2002.67E-05
ENSG00000126353TMCO10.647-0.435-3.9697.22E-05
ENSG00000115085MAP2K70.649-0.432-3.9169.01E-05
ENSG00000170482SPINK10.538-0.620-3.8969.80E-05
ENSG00000162545GLYCTK1.5030.4073.8800.0001
ENSG00000125910ERRFI10.638-0.449-3.8380.0001
ENSG00000174652MSANTD30.674-0.395-3.7950.0001
ENSG00000197540BTLA0.655-0.423-3.7220.0001
ENSG00000184903ATP6V0E11.3920.3313.7190.0001
ENSG00000132465CALML50.692-0.369-3.6980.0002
ENSG00000089199MORF4L21.3820.3243.6830.0002
ENSG00000189319ZZEF10.696-0.362-3.6740.0002
ENSG00000198198TYK20.692-0.369-3.6090.0003
ENSG00000153531EOMES1.3180.2763.5530.0003
ENSG00000159176PTGR11.3590.3073.5520.0003
ENSG00000127241GPR1710.598-0.515-3.5370.0004
ENSG00000197992LIMD20.617-0.483-3.5370.0004
ENSG00000182866DCP20.700-0.357-3.4940.0004

Recognition of gene sets for genome variation

Using copy number variation data in TCGA, we used GISTIC 2.0 to identify genes exhibiting significant amplification or deletion. Significantly amplified fragments within the genomes, including EGFR at 7p11.2 (q value = 2.28E-43), FGF1 at 8p11.23 (q value = 3.94E-14), and ERBB2 at 17q12 (q value = 0.0050541) (Figure 1A). In total, 247 genes were amplified (Supplementary Table 1). Genes that were significantly missing from the genome included CDKN2A at 9p21.3 (q value = 5.28E-149), CDK5 at 7q36.1 (q value = 1.87E-06), and PTEN at 10q23.31 (q = 0.0032849) (Figure 1B). A total of 901 genes were identified as missing from the genome.

Identification of genes with significant amplification or deletion. (A and B) The mRNA located in the focal CNA peaks are HNSCC-related. False-discovery rates (q values) and scores from GISTIC 2.0 for alterations (x-axis) are plotted against genome positions (y-axis). Dotted lines indicate the centromeres. Amplifications (A) are shown in red, deletion (B) in blue. The green line represents 0.25 q value cut-off point that determines significance. (C) Top 50 genes with the most significant mutations. The bar chart above shows the total number of synonymous and non-synonymous mutations in each patient's top 50 genes. The bar chart on the right shows the number of samples in which the 50 genes were mutated in all samples. The different colors in the thermogram indicate the type of mutation; gray indicates no mutation.

Figure 1. Identification of genes with significant amplification or deletion. (A and B) The mRNA located in the focal CNA peaks are HNSCC-related. False-discovery rates (q values) and scores from GISTIC 2.0 for alterations (x-axis) are plotted against genome positions (y-axis). Dotted lines indicate the centromeres. Amplifications (A) are shown in red, deletion (B) in blue. The green line represents 0.25 q value cut-off point that determines significance. (C) Top 50 genes with the most significant mutations. The bar chart above shows the total number of synonymous and non-synonymous mutations in each patient's top 50 genes. The bar chart on the right shows the number of samples in which the 50 genes were mutated in all samples. The different colors in the thermogram indicate the type of mutation; gray indicates no mutation.

Using TCGA mutation annotation data with Mutsig2, we identified genes with significant mutations. A total of 302 genes with significant mutation frequencies were detected (Supplementary Table 2). The most significant types mutations in the Top 50 genes were synonymous mutations, missense mutations, frame insertions or deletions, frame movement, nonsense mutations, distribution of shear sites, and other non-synonymous mutations (Figure 1C). It was clear that there were differences in the mutations to different genes, including B2M, CDKN2A, PTEN, TP53 and FBXW7. A common feature, however, is that mutation of these genes is reported to be closely related to the occurrence and development of tumors [7–10].

Functional analysis of genome variant genes

To analyze the function of genomic variant genes, we integrated 1321 amplified or deleted genes and significantly mutated genes identified based on copy number variation. Gene Ontology (GO) biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were then performed. The results of the KEGG enrichment analysis revealed that Pathways in cancer, HPV infection, PI3K-Akt signaling pathway, Human T-cell leukemia virus 1 infection, Human cytomegalovirus infection, and numerous other related KEGG biological pathways are important to the development of cancer (Figure 2A). GO terms such as developmental process, positive regulation of cellular process, cell differentiation, and regulation of localization were mainly enriched in the “biological process” category (Figure 2B). These terms were also closely related to the occurrence and development of cancer; that is, these genes exhibiting genomic variation are closely related to cancer.

Functional enrichment analysis of 1321 genome variant genes. (A) Enriched KEGG biological pathways. (B) Enriched GO terms in the “biological process” category. Different colors indicate different significances, while different sizes indicate the number of genes.

Figure 2. Functional enrichment analysis of 1321 genome variant genes. (A) Enriched KEGG biological pathways. (B) Enriched GO terms in the “biological process” category. Different colors indicate different significances, while different sizes indicate the number of genes.

Identification of a 6-gene signature for head and neck cancer survival

To identify a gene signature, we integrated genomic variant and prognosis-related genes, and then selected the intersection of the groups as candidate genes, which yielded 36 genes. We then used random forests for feature selection. The relationship between error rate and number of taxonomic trees was used to reveal genes with relative importance greater than 0.4 as the final signature (Figure 3A). Ultimately, we identified 6 genes (Table 2). The important order of the out-of-bag scores for the 6 genes is displayed in Figure 3B. A 6-gene signature was established by multivariate COX regression analysis. The Equation 1 is as follows:

Riska =-0.3929517* D2HGDH+ 0.3627132 * PEX11A + 0.3038125 * NLRP2 + 0.275704 * SERP1NE1 + 0.188539* UPK2+0.1112888* CTTN

Identification of genomic variant genes and prognosis-related genes in head and neck cancer. (A) Relationship between the error rate and number of classification trees. (B) Importance the sequencing of 6 out-of-bag genes. (C) Distribution of the 6-gene signature in Kaplan-Meier survival curve for the TCGA training set. (D) ROC curve and AUC for the 6-gene signature classification. (E) TCGA training focused on risk score, survival time, survival status, and expression of the 6-gene signature.

Figure 3. Identification of genomic variant genes and prognosis-related genes in head and neck cancer. (A) Relationship between the error rate and number of classification trees. (B) Importance the sequencing of 6 out-of-bag genes. (C) Distribution of the 6-gene signature in Kaplan-Meier survival curve for the TCGA training set. (D) ROC curve and AUC for the 6-gene signature classification. (E) TCGA training focused on risk score, survival time, survival status, and expression of the 6-gene signature.

Table 2. Six genes significantly associated with OS in the training set patients.

Ensemble Gene IDSymbolHRZ-scorePImportanceRelative importance
ENSG00000180902D2HGDH0.77-2.618.88E-030.0121
ENSG00000166821PEX11A1.332.982.87E-030.00760.6359
ENSG00000022556NLRP21.322.963.03E-030.00740.6214
ENSG00000106366SERPINE11.433.338.47E-040.00520.4369
ENSG00000110375UPK21.292.982.84E-030.0050.4175
ENSG00000085733CTTN1.272.697.08E-030.00480.4

The risk score of each sample was calculated, and the samples were grouped according to the median risk score (cutoff = -0.0236503). The prognoses of the high-risk and low-risk groups significantly differed (Figure 3C). The average 1-, 3-, and 5-year AUC for the 6-gene signature was 0.75 (Figure 3D). High expression of PEX11A, NLRP2, SERPINE1, UPK2 and CTTN was associated with high risk, while high expression of D2HGDH was associated with low risk and was a protective factor.

Verification of the robustness of the 6-gene signature model

To verify the robustness of the 6-gene signature model, we first calculated a risk score for each sample in the test set. Based on the threshold for the training set, the samples were divided into two groups with significantly different prognoses (Figure 4A). The relationship between the expression of the 6 genes and the risk score was also consistent with the training set (Figure 4C). We similarly applied the model to all TCGA tumor samples and found that the low risk group fared significantly better than the high-risk group (Figure 4B). The relationship between the expression of these 6 genes and the risk score was also consistent with the training set (Figure 4D). Overall, the model effectively provided prognostic classifications with the TCGA datasets.

Relation between the 6-gene signature and cancer risk. (A) Kaplan-Meier curve for the test set sample. (B) Kaplan-Meier curve in all TCGA tumor samples. (C) Relationship between expression of the 6-gene signature and risk scores in test set samples. (D) Relationship between expression of the 6-gene signature and the risk score in all TCGA samples.

Figure 4. Relation between the 6-gene signature and cancer risk. (A) Kaplan-Meier curve for the test set sample. (B) Kaplan-Meier curve in all TCGA tumor samples. (C) Relationship between expression of the 6-gene signature and risk scores in test set samples. (D) Relationship between expression of the 6-gene signature and the risk score in all TCGA samples.

To verify the classification performance of the 6-gene signature model with different data platforms, we used GEO platform data as external datasets, used the model to calculate a risk score for each sample, and used the cutoff for the training set to divide the samples into high-risk and low-risk groups. The prognosis of the low-risk group was significantly better than that of the high-risk group (Figure 5A). ROC analysis showed that the 5-year AUC was up to 0.74, compared with the training set. As in Figure 5B, the data in Figure 5C show that the relationship between the expression of the 6 genes and risk score is also consistent with the training set. Thus, the 6-gene signature model we selected was prognostic with both internal and external datasets.

Performance of the 6-gene signature model with GEO data. (A) Kaplan-Meier survival curve distribution of 6-gene signature for the GSE65858 dataset. (B) ROC curve and AUC for the 6-gene signature classification. (C) Risk score, survival time, survival status, and expression of the 6-gene signature in the GSE65858 dataset.

Figure 5. Performance of the 6-gene signature model with GEO data. (A) Kaplan-Meier survival curve distribution of 6-gene signature for the GSE65858 dataset. (B) ROC curve and AUC for the 6-gene signature classification. (C) Risk score, survival time, survival status, and expression of the 6-gene signature in the GSE65858 dataset.

Clinical independence of the 6-gene signature model

To assess the independence of the 6-gene signature model in clinical application, we used single-factor and multi-factor COX regression to analyze HRs, 95% CIs, and P values from TCGA training set, TCGA test set, and the GSE65858 data. We systematically analyzed the clinical information from the patients as recorded in TCGA and GSE65858 datasets, including their age, sex, disease stage, pathological TNM stage, and tumor stage, as well as our 6-gene signature (Table 3). In TCGA test set, single factor COX regression analysis revealed that in the high-risk group, pathologic T3, T4, M1, and N2 all significantly correlated with survival. However, the corresponding multifactor COX regression analysis found that only high-risk group (HR = 2.17, 95% CI = 1.04-4.53, P = 0.037979), pathologic T4 (HR = 7.58, 95% CI = 1.90-30.22, P = 0.004054) and pathologic N2 (HR = 5.01, 95% CI = 1.97-12.75, P = 0.000703) had clinical independence. In TCGA test set, single factor COX regression analysis revealed that in the high-risk group, pathologic T3 and N2 correlated significantly with survival, but the corresponding multivariate COX regression analysis found that no factor had clinical independence; the HR in the high-risk group was 1.45, 95% CI = 0.67-3.12, P = 0.3389. In the GSE65858 dataset, univariate COX regression analysis revealed that the in the high-risk group, age, pathologic T3, N3, and M1 correlated significantly with survival. The corresponding multivariate COX regression analysis found that the high-risk group (HR = 1.83, 95% CI = 1.16-2.87, P = 0.0083) and age (HR = 1.03, 95% CI = 1.01-1.05, P = 0.00313) had clinical independence. These results show that our model 6-gene signature is a prognostic index independent of other clinical factors and exhibits independent predictive performance upon clinical application.

Table 3. Univariate and multivariate COX regression analyses of clinical factors and independence associated with prognosis.

VariablesUnivariate analysisMultivariable analysis
HR95%CI of HRPHR95%CI of HRP
TCGA training datasets
6-gene risk score
Low risk group1(reference)1(reference)
High risk group2.661.79-3.921.060E-062.181.04-4.530.038
Age1.031.01-1.048.550E-040.990.95-1.020.604231
Gender female1(reference)1(reference)
Gender male0.730.49- 1.070.110.500.23-1.090.081554
Grade 11(reference)1(reference)
Grade 21.780.95-3.320.070.910.25-3.180.88025
Grade 3 / 41.470.75-2.870.261.940.53-7.040.315185
Pathologic T 1/ T 21(reference)1(reference)
Pathologic T 32.261.28-3.954.53E-032.070.51-8.370.309
Pathologic T 42.461.50-4.003.19E-047.591.90-30.220.004
Pathologic N 01(reference)1(reference)
Pathologic N 10.890.37-2.080.7820.760.19-3.010.694
Pathologic N 22.291.40-3.720.0015.021.97-12.750.001
Pathologic M 01(reference)1(reference)
Pathologic M 1/ M X1.310.66-2.550.4331.160.43-3.080.76
Tumor stage I1(reference)1(reference)
Tumor stage II1.080.23-4.900.9180.550.03-9.850.686335
Tumor stage III1.400.30-6.420.6670.770.05-10.600.845959
Tumor stage IV2.750.67-11.220.1580.270.02-3.630.32
Validation cohort, TCGA test datasets and GSE65858
TCGA test datasets
6-gene risk score
Low risk group1(reference)1(reference)
High risk group1.621.08- 2.420.0201.450.67-3.120.339
Age1.010.98-1.020.4441.000.96-1.030.993
Gender female1(reference)1(reference)
Gender male0.840.55-1.280.4200.450.19-1.010.054
Grade 11(reference)1(reference)
Grade 21.860.92-3.770.080.900.28-2.880.865
Grade 31.570.73-3.360.241.540.42-5.520.508
Pathologic T 1/ T 21(reference)1(reference)
Pathologic T 31.841.09-3.110.0221.760.47-6.540.399
Pathologic T 41.360.84-2.210.2131.200.35-4.070.770
Pathologic N 01(reference)
Pathologic N 11.140.57-2.250.7061.740.50-60.379
Pathologic N 22.491.50-4.120.0001.870.65-5.320.239
Pathologic N 32.900.87-9.60.0821.860.26-12.820.529
Pathologic M 01(reference)1(reference)
Pathologic M 11.050.49-2.210.9051.340.54-3.280.524
Tumor stage I1(reference)
Tumor stage II2.940.34-0.660.1570.550.03-9.200.677
Tumor stage III3.040.32-0.700.1371.000.08-11.420.997
Tumor stage IV4.030.24-0.980.0531.460.12-17.690.765
GSE65858
6-gene risk score
Low risk group1(reference)1(reference)
High risk group1.751.13-2.670.0111.831.16-2.870.008
Age1.031.00-1.040.011.031.01-1.050.00
Gender female1(reference)1(reference)
Gender male1.050.61-1.770.871.020.59-1.750.94
Pathologic T 11(reference)1(reference)
Pathologic T 20.490.21-1.150.100.530.15-1.780.31
Pathologic T 31.730.81-3.670.151.610.55-4.690.38
Pathologic T 41.890.91-3.870.081.220.42-3.530.71
Pathologic N 01(reference)1(reference)
Pathologic N 10.360.12-1.020.060.340.10-1.080.07
Pathologic N 21.600.99-2.560.051.030.51-2.050.94
Pathologic N 32.931.34-6.420.011.270.41-3.840.67
Pathologic M 01(reference)1(reference)
Pathologic M 1/ M X3.711.63-8.40.002.180.70-6.790.18
Tumor stage I1(reference)1(reference)
Tumor stage II0.390.11-1.330.130.600.10-3.400.56
Tumor stage III0.460.14-1.450.190.660.13-3.310.61
Tumor stage IV1.490.60-3.700.391.140.24-5.200.87

GSEA analysis of pathway differences enriched in the high-risk and low-risk groups

GSEA was used in TCGA training to analyze the pathways significantly enriched in the high-risk and low-risk groups. Twenty enriched pathways were detected (Supplementary Table 3), including focal adhesion, TGF-β signaling pathway, WNT signaling pathway, and ERBB signaling pathway, all of which are closely related to tumor occurrence, development and metastasis. Notably, these pathways are significantly enriched in the high-risk samples (Figure 6).

GSEA showing four pathways enriched in the high-risk group. GSEA enrichment results for focal adhesion, TGF-β signaling pathway, WNT signaling pathway, and ERBB signaling pathway.

Figure 6. GSEA showing four pathways enriched in the high-risk group. GSEA enrichment results for focal adhesion, TGF-β signaling pathway, WNT signaling pathway, and ERBB signaling pathway.

Discussion

In terms of prognosis, head and neck cancer is a highly heterogeneous disease in that survival times vary substantially among patients with similar TNM stages. With the diagnosis and treatment of head and neck cancers at earlier stages, traditional clinicopathological indicators such as tumor size, vascular invasion, portal vein thrombus and TNM stage have proven inadequate for predicting individual outcomes, especially risk stratification, as no one-size-fits-all treatment strategy appears to be effective [11, 12]. Consequently, screening prognostic molecular markers that adequately reflect the biological characteristics of tumors would be of great significance for individualized prevention and treatment of head and neck cancer patients. In the present study, we analyzed the expression profiles of 771 head and neck cancer samples from TCGA and the GEO and identified 6 genes robustly associated with OS. This signature is independent of other clinical factors.

Gene signatures are currently being used in clinical practice. Two examples are Oncotype DX [13–15], which provides a breast cancer recurrence score based on expression of 21 genes, and Coloprint, which provides a colon cancer recurrence score based on expression of 18 genes [16–18]. Results obtained with these assays have shown that screening new prognostic cancer markers based on gene expression profiles is a promising high-throughput molecular identification method. In that regard, Tian et al. [6] identified a 6-gene signature, but the verification set AUC was only about 0.65, and Zhao et al. [4] identified a 3-gene signature, but the AUC was only about 0.6. In addition, De et al. 5] identified a 172-gene signature through meta-analysis of gene expression. Although the AUC is high, the large number of genes that needs to be detected makes this analysis impractical for clinical use. By contrast, our 6-gene signature has a high AUC using only 6 genes, which makes it conducive to clinical application.

The six genes in our signature include PEX11A, NLRP2, SERPINE1, UPK, and CTTN as risk factors, and D2HGDH as a protective factor. It has been reported that NLRP2 can be used as a marker of leukemia [19] and has an important impact on the prognosis after stem cell transplantation [20]. SERPINE1 is closely related to prognosis in ovarian cancer, gastric cancer, thyroid cancer and other tumors [21–25]. CTTN is closely related to prognosis in head and neck cancer, esophageal cancer, thyroid cancer, and glioma, among others [26–29]. D2HGDH is a marker of colorectal cancer, glioma, and prostate cancer [30–33]. PEX11A and UPK have not been previously reported to be related to cancer. Ours is the first study to suggest that they can be used as new prognostic markers of head and neck cancer. At the same time, our GSEA results show that the 6-gene signature enrichment significantly correlates with pathways and biological processes associated with the occurrence and development of HNSCC. This indicates that our model has potential clinical application value and could provide a potential target for diagnosis and for development of new targeted therapies that include, for example, novel alkylating agents [34, 35].

Although we have identified potential candidate genes affecting tumor prognosis using bioinformatics technology with large samples, our study has limitations. First, the sample lacks some clinical follow-up information, so we did not consider factors such as the presence of other health conditions to differentiate prognostic biomarkers. Second, the results obtained using bioinformatics analysis alone are insufficient and need to be confirmed through experimental verification. Therefore, further genetic and experimental studies with larger samples and experimental validation are needed.

Conclusions

In our study, we developed a 6-gene signature prognostic stratification system, which has good AUC values in both the training set and validation set, and is independent of other clinical features. Compared to clinical features, gene classifiers can improve survival risk prediction. We therefore recommend using this classifier as a molecular diagnostic test to assess prognostic risk in patients with head and neck cancer.

Materials and Methods

Data acquisition and processing

The FPKM data of TCGA RNA-Seq downloaded from the UCSC Cancer Browser (https://xenabrowser.net/datapages/) contained 546 samples. The clinical follow-up information contained 612 samples, and the copy number variation data on the SNP 6.0 chip contained 519 samples. The mutation annotation information (MAF) downloaded from the GDC client contained 504 samples. The standardized tables of the GSE658587 [36] data set were downloaded from the GEO. For TCGA RNAseq data, 501 tumor samples with follow-up information were chosen and randomly divided into two groups: a training set (N = 250) (Supplementary Tables 4 and 5) and a test set (N = 251) (Supplementary Tables 6 and 7). The GSE65858 data set was used as an external verification set for each group of sample information (Table 4).

Table 4. Clinical information statistics for three datasets.

CharacteristicTCGA training datasets (n=250)TCGA test datasets (n=251)GSE65858 (n=270)
Age(years)<=50385041
>50212201229
Survival StatusLiving136146176
Dead14110594
Genderfemale686647
male182185223
GradeG 13428
G 2136163
G 36851
G 420
Pathologic_TT 1172835
T 2676580
T 3484858
T 4896397
Pathologic_NN 0809094
N 1263932
N 29472132
N 32512
Pathologic_MM 09493263
M 1/ M X32307
Tumor StageStage I91618
Stage II373237
Stage III304837
Stage IV141120178

Univariate Cox proportional risk regression analysis

As described previously by Jin-Cheng et al. [37], univariate Cox proportional hazard regression analysis was performed with each immune gene to screen out genes significantly associated with OS in the training data set. P < 0.01 was chosen as the threshold.

Analysis of copy number variation data

GISTIC software is widely used to detect both broad and focal (potentially overlapping) recurring events. We used GISTIC 2.0 [38] to identify genes exhibiting significant amplification or deletion. The parameter threshold was that the length of the amplification or deletion was more than 0.1, and the fragment was P < 0.05.

Gene mutation analysis

We used Mutsig 2.0 software with TCGA mutation data to identify genes with significant mutations in their MAF files. The threshold used was P < 0.05.

Construction of a prognostic gene signature

We chose the genes that were significantly related to OS and genes that were exhibited amplification, deletion or mutation. We further used the Random Survival Forest algorithm to rank the importance of prognostic genes. Like Jin et al. [39], we used the R package random Survival Forest to screen the prognostic genes. We set the number of Monte Carlo iterations to 100 and the number of steps forward to 5, and identified the genes whose relative importance as characteristic genes was greater than 0.4. In addition, we carried out a multivariate Cox regression analysis, and the following risk scoring model was constructed using the Equation 2:

RickScore​ =∑k−lnExpk∗eHRk

where N is the number of prognostic genes, Expk is the expression value of the prognostic genes, and eHRk is the estimated regression coefficient of genes in the multivariate Cox regression analysis.

Functional enrichment analyses

GO and KEGG pathway enrichment analysis was performed using the R package cluster profiler for genes [40], to identify over-represented GO terms in three categories (biological processes, molecular function and cellular component) as well as over-represented KEGG pathway terms. For this analyses, a false discovery rate < 0.05 was considered to indicate statistical significance.

GSEA [41] was performed using the http://software.broadinstitute.org/gsea/downloads.jsp website with the MSigDB [42] C2 Canonical pathways gene set collection, which contains 1320 gene sets. Gene sets with a false discovery rate < 0.05 after performing 1000 permutations were considered to be significantly enriched.

Statistical analysis

When the median risk score in each data set was used as a cutoff to compare survival risk between high-risk and low-risk groups, a Kaplan-Meier (KM) curve was drawn. Multivariate Cox regression analysis was used to test whether gene markers were independent prognostic factors. Significance was defined as P < 0.05. All analyses were performed using R 3.4.3.

Abbreviations

HNSCC: head and neck squamous cell carcinoma; TCGA: the cancer genome atlas; GEO: group on earth observations; SNPs: single nucleotide polymorphisms; CNV: copy number variation; SCC-A: squamous cell carcinoma antigen; HPV: human papilloma virus; MFA: the mutation annotation information; OS: overall survival; GO: gene ontology; KEGG: kyoto encyclopedia of genes and genomes; KM: Kaplan-Meier; HR: higher risk; PPI: protein-protein interaction.

Conflicts of Interest

The authors declare that they have no competing interests.

Funding

The work was supported by The Natural Science Foundation of China (no.81660294, 81602389, 81472696, 81202128, 81272974) and The Natural Science Youth Founds of Jiangxi Province (no.2016BAB215249).

References

  • 1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2016. CA Cancer J Clin. 2016; 66:7–30. https://doi.org/10.3322/caac.21332 [PubMed]
  • 2. Leemans CR, Braakhuis BJ, Brakenhoff RH. The molecular biology of head and neck cancer. Nat Rev Cancer. 2011; 11:9–22. https://doi.org/10.1038/nrc2982 [PubMed]
  • 3. Chauhan SS, Kaur J, Kumar M, Matta A, Srivastava G, Alyass A, Assi J, Leong I, MacMillan C, Witterick I, Colgan TJ, Shukla NK, Thakar A, et al. Prediction of recurrence-free survival using a protein expression-based risk classifier for head and neck cancer. Oncogenesis. 2015; 4:e147. https://doi.org/10.1038/oncsis.2015.7 [PubMed]
  • 4. Zhao X, Sun S, Zeng X, Cui L. Expression profiles analysis identifies a novel three-mRNA signature to predict overall survival in oral squamous cell carcinoma. Am J Cancer Res. 2018; 8:450–61. [PubMed]
  • 5. De Cecco L, Bossi P, Locati L, Canevari S, Licitra L. Comprehensive gene expression meta-analysis of head and neck squamous cell carcinoma microarray data defines a robust survival predictor. Ann Oncol. 2014; 25:1628–35. https://doi.org/10.1093/annonc/mdu173 [PubMed]
  • 6. Tian S, Meng G, Zhang W. A six-mRNA prognostic model to predict survival in head and neck squamous cell carcinoma. Cancer Manag Res. 2018; 11:131–42. https://doi.org/10.2147/CMAR.S185875 [PubMed]
  • 7. Yeon Yeon S, Jung SH, Jo YS, Choi EJ, Kim MS, Chung YJ, Lee SH. Immune checkpoint blockade resistance-related B2M hotspot mutations in microsatellite-unstable colorectal carcinoma. Pathol Res Pract. 2019; 215:209–14. https://doi.org/10.1016/j.prp.2018.11.014 [PubMed]
  • 8. Padhi SS, Roy S, Kar M, Saha A, Roy S, Adhya A, Baisakh M, Banerjee B. Role of CDKN2A/p16 expression in the prognostication of oral squamous cell carcinoma. Oral Oncol. 2017; 73:27–35. https://doi.org/10.1016/j.oraloncology.2017.07.030 [PubMed]
  • 9. Hou SQ, Ouyang M, Brandmaier A, Hao H, Shen WH. PTEN in the maintenance of genome integrity: from DNA replication to chromosome segregation. BioEssays. 2017; 39:1700082. https://doi.org/10.1002/bies.201700082 [PubMed]
  • 10. Soussi T, Wiman KG. TP53: an oncogene in disguise. Cell Death Differ. 2015; 22:1239–49. https://doi.org/10.1038/cdd.2015.53 [PubMed]
  • 11. Llovet JM, Ricci S, Mazzaferro V, Hilgard P, Gane E, Blanc JF, de Oliveira AC, Santoro A, Raoul JL, Forner A, Schwartz M, Porta C, Zeuzem S, et al, and SHARP Investigators Study Group. Sorafenib in advanced hepatocellular carcinoma. N Engl J Med. 2008; 359:378–90. https://doi.org/10.1056/NEJMoa0708857 [PubMed]
  • 12. Cheng AL, Kang YK, Chen Z, Tsao CJ, Qin S, Kim JS, Luo R, Feng J, Ye S, Yang TS, Xu J, Sun Y, Liang H, et al. Efficacy and safety of sorafenib in patients in the Asia-Pacific region with advanced hepatocellular carcinoma: a phase III randomised, double-blind, placebo-controlled trial. Lancet Oncol. 2009; 10:25–34. https://doi.org/10.1016/S1470-2045(08)70285-7 [PubMed]
  • 13. Siow ZR, De Boer RH, Lindeman GJ, Mann GB. Spotlight on the utility of the Oncotype DX® breast cancer assay. Int J Womens Health. 2018; 10:89–100. https://doi.org/10.2147/IJWH.S124520 [PubMed]
  • 14. Bhutiani N, Egger ME, Ajkay N, Scoggins CR, Martin RC 2nd, McMasters KM. Multigene signature panels and breast cancer therapy: patterns of use and impact on clinical decision making. J Am Coll Surg. 2018; 226:406–412.e1. https://doi.org/10.1016/j.jamcollsurg.2017.12.043 [PubMed]
  • 15. Wang SY, Dang W, Richman I, Mougalian SS, Evans SB, Gross CP. cost-effectiveness analyses of the 21-gene assay in breast cancer: systematic review and critical appraisal. J Clin Oncol. 2018; 36:1619–27. https://doi.org/10.1200/JCO.2017.76.5941 [PubMed]
  • 16. Kopetz S, Tabernero J, Rosenberg R, Jiang ZQ, Moreno V, Bachleitner-Hofmann T, Lanza G, Stork-Sloots L, Maru D, Simon I, Capellà G, Salazar R. Genomic classifier ColoPrint predicts recurrence in stage II colorectal cancer patients more accurately than clinical factors. Oncologist. 2015; 20:127–33. https://doi.org/10.1634/theoncologist.2014-0325 [PubMed]
  • 17. Tan IB, Tan P. Genetics: an 18-gene signature (ColoPrint®) for colon cancer prognosis. Nat Rev Clin Oncol. 2011; 8:131–33. https://doi.org/10.1038/nrclinonc.2010.229 [PubMed]
  • 18. Maak M, Simon I, Nitsche U, Roepman P, Snel M, Glas AM, Schuster T, Keller G, Zeestraten E, Goossens I, Janssen KP, Friess H, Rosenberg R. Independent validation of a prognostic genomic signature (ColoPrint) for patients with stage II colon cancer. Ann Surg. 2013; 257:1053–58. https://doi.org/10.1097/SLA.0b013e31827c1180 [PubMed]
  • 19. Zhao X, Li Y, Wu H. A novel scoring system for acute myeloid leukemia risk assessment based on the expression levels of six genes. Int J Mol Med. 2018; 42:1495–507. https://doi.org/10.3892/ijmm.2018.3739 [PubMed]
  • 20. Granell M, Urbano-Ispizua A, Pons A, Aróstegui JI, Gel B, Navarro A, Jansa S, Artells R, Gaya A, Talarn C, Fernández-Avilés F, Martínez C, Rovira M, et al. Common variants in NLRP2 and NLRP3 genes are strong prognostic factors for the outcome of HLA-identical sibling allogeneic stem cell transplantation. Blood. 2008; 112:4337–42. https://doi.org/10.1182/blood-2007-12-129247 [PubMed]
  • 21. Pan JX, Qu F, Wang FF, Xu J, Mu LS, Ye LY, Li JJ. Aberrant SERPINE1 DNA methylation is involved in carboplatin induced epithelial-mesenchymal transition in epithelial ovarian cancer. Arch Gynecol Obstet. 2017; 296:1145–52. https://doi.org/10.1007/s00404-017-4547-x [PubMed]
  • 22. Azimi I, Petersen RM, Thompson EW, Roberts-Thomson SJ, Monteith GR. Hypoxia-induced reactive oxygen species mediate N-cadherin and SERPINE1 expression, EGFR signalling and motility in MDA-MB-468 breast cancer cells. Sci Rep. 2017; 7:15140. https://doi.org/10.1038/s41598-017-15474-7 [PubMed]
  • 23. Yu XM, Jaskula-Sztul R, Georgen MR, Aburjania Z, Somnay YR, Leverson G, Sippel RS, Lloyd RV, Johnson BP, Chen H. Notch1 signaling regulates the aggressiveness of differentiated thyroid cancer and inhibits SERPINE1 expression. Clin Cancer Res. 2016; 22:3582–92. https://doi.org/10.1158/1078-0432.CCR-15-1749 [PubMed]
  • 24. Trabert B, Eldridge RC, Pfeiffer RM, Shiels MS, Kemp TJ, Guillemette C, Hartge P, Sherman ME, Brinton LA, Black A, Chaturvedi AK, Hildesheim A, Berndt SI, et al. Prediagnostic circulating inflammation markers and endometrial cancer risk in the prostate, lung, colorectal and ovarian cancer (PLCO) screening trial. Int J Cancer. 2017; 140:600–10. https://doi.org/10.1002/ijc.30478 [PubMed]
  • 25. Mazzoccoli G, Pazienza V, Panza A, Valvano MR, Benegiamo G, Vinciguerra M, Andriulli A, Piepoli A. ARNTL2 and SERPINE1: potential biomarkers for tumor aggressiveness in colorectal cancer. J Cancer Res Clin Oncol. 2012; 138:501–11. https://doi.org/10.1007/s00432-011-1126-6 [PubMed]
  • 26. Ramos-García P, González-Moles MA, Ayén Á, González-Ruiz L, Ruiz-Ávila I, Gil-Montoya JA. Prognostic and clinicopathological significance of CTTN/cortactin alterations in head and neck squamous cell carcinoma: systematic review and meta-analysis. Head Neck. 2019; 41:1963–78. https://doi.org/10.1002/hed.25632 [PubMed]
  • 27. Zhang S, Qi Q. MTSS1 suppresses cell migration and invasion by targeting CTTN in glioblastoma. J Neurooncol. 2015; 121:425–31. https://doi.org/10.1007/s11060-014-1656-2 [PubMed]
  • 28. Hu X, Moon JW, Li S, Xu W, Wang X, Liu Y, Lee JY. Amplification and overexpression of CTTN and CCND1 at chromosome 11q13 in Esophagus squamous cell carcinoma (ESCC) of North Eastern Chinese Population. Int J Med Sci. 2016; 13:868–74. https://doi.org/10.7150/ijms.16845 [PubMed]
  • 29. Long HC, Gao X, Lei CJ, Zhu B, Li L, Zeng C, Huang JB, Feng JR. miR-542-3p inhibits the growth and invasion of colorectal cancer cells through targeted regulation of cortactin. Int J Mol Med. 2016; 37:1112–18. https://doi.org/10.3892/ijmm.2016.2505 [PubMed]
  • 30. Han J, Jackson D, Holm J, Turner K, Ashcraft P, Wang X, Cook B, Arning E, Genta RM, Venuprasad K, Souza RF, Sweetman L, Theiss AL. Elevated d-2-hydroxyglutarate during colitis drives progression to colorectal cancer. Proc Natl Acad Sci USA. 2018; 115:1057–62. https://doi.org/10.1073/pnas.1712625115 [PubMed]
  • 31. Qin F, Song Z, Chang M, Song Y, Frierson H, Li H. Recurrent cis-SAGe chimeric RNA, D2HGDH-GAL3ST2, in prostate cancer. Cancer Lett. 2016; 380:39–46. https://doi.org/10.1016/j.canlet.2016.06.013 [PubMed]
  • 32. Brehmer S, Pusch S, Schmieder K, von Deimling A, Hartmann C. Mutational analysis of D2HGDH and L2HGDH in brain tumours without IDH1 or IDH2 mutations. Neuropathol Appl Neurobiol. 2011; 37:330–32. https://doi.org/10.1111/j.1365-2990.2010.01114.x [PubMed]
  • 33. Krell D, Assoku M, Galloway M, Mulholland P, Tomlinson I, Bardella C. Screen for IDH1, IDH2, IDH3, D2HGDH and L2HGDH mutations in glioblastoma. PLoS One. 2011; 6:e19868. https://doi.org/10.1371/journal.pone.0019868 [PubMed]
  • 34. Kou Y, Koag MC, Lee S. N7 methylation alters hydrogen-bonding patterns of guanine in duplex DNA. J Am Chem Soc. 2015; 137:14067–70. https://doi.org/10.1021/jacs.5b10172 [PubMed]
  • 35. Kou Y, Koag MC, Lee S. Structural and kinetic studies of the effect of guanine N7 alkylation and metal cofactors on DNA replication. Biochemistry. 2018; 57:5105–16. https://doi.org/10.1021/acs.biochem.8b00331 [PubMed]
  • 36. Wichmann G, Rosolowski M, Krohn K, Kreuz M, Boehm A, Reiche A, Scharrer U, Halama D, Bertolini J, Bauer U, Holzinger D, Pawlita M, Hess J, et al, and Leipzig Head and Neck Group (LHNG). The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer. Int J Cancer. 2015; 137:2846–57. https://doi.org/10.1002/ijc.29649 [PubMed]
  • 37. Guo JC, Wu Y, Chen Y, Pan F, Wu ZY, Zhang JS, Wu JY, Xu XE, Zhao JM, Li EM, Zhao Y, Xu LY. Protein-coding genes combined with long noncoding RNA as a novel transcriptome molecular staging model to predict the survival of patients with esophageal squamous cell carcinoma. Cancer Commun (Lond). 2018; 38:4. https://doi.org/10.1186/s40880-018-0277-0 [PubMed]
  • 38. Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011; 12:R41. https://doi.org/10.1186/gb-2011-12-4-r41 [PubMed]
  • 39. Meng J, Li P, Zhang Q, Yang Z, Fu S. A four-long non-coding RNA signature in predicting breast cancer survival. J Exp Clin Cancer Res. 2014; 33:84. https://doi.org/10.1186/s13046-014-0084-7 [PubMed]
  • 40. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
  • 41. Subramanian A, Kuehn H, Gould J, Tamayo P, Mesirov JP. GSEA-P: a desktop application for Gene Set Enrichment Analysis. Bioinformatics. 2007; 23:3251–53. https://doi.org/10.1093/bioinformatics/btm369 [PubMed]
  • 42. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011; 27:1739–40. https://doi.org/10.1093/bioinformatics/btr260 [PubMed]
Download PDF
Corresponding Authors

Gengming Cai
cgmkgx@hotmail.com

Yunxia Lv
83394045@qq.com

Keywords
bioinformatics CNV prognostic markers TCGA
Paper Sections