Non-negative matrix factorization model-based construction for molecular clustering and prognostic assessment of head and neck squamous carcinoma.

Purpose: We aimed at exploring the efficacy of non-negative matrix factorization (NMF) model-based clustering for prognostic assessment of head and neck squamous carcinoma (HNSCC).
Methods: The transcriptome microarray data of HNSCC samples were downloaded from The Cancer Genome Atlas (TCGA) and the Shanghai Ninth People's Hospital. R software packages were used to establish NMF clustering, from which relevant prognostic models were developed.
Results: Based on NMF, samples were allocated into 2 subgroups. Predictive models were constructed using differentially expressed genes between the two subgroups. The high-risk group was associated with poor prognostic outcomes. Moreover, multi-factor Cox regression analysis revealed that the predictive model was an independent prognostic predictor.
Conclusion: The NMF-based prognostic model has the potential for prognostic assessment of HNSCC.

Introduction

Globally, head and neck squamous cell carcinoma (HNSCC) is one of the most common malignant tumors, accounting for more than 90% of all malignant tumors of the head and neck [1]. Treatment options for HNSCC are mainly based on TNM staging and a combination of surgical-based therapies (radiotherapy, chemotherapy and biotherapy) [2]. Although a majority of HNSCC patients present with locally advanced disease with significant lymph node metastases, advances in multi-disciplinary treatment has improved treatment outcomes. However, HNSCC-associated mortality rate is still above 55%, with 40–60% recurrence and metastasis rates [2, 3, 4]. Therefore, accurate prognostic prediction of HNSCC patients is important for informing clinical treatment. Clustering of HNSCC samples and analysis of biological differences between groups are useful for elucidating the pathogenesis of HNSCC. They also have a reference value for clinical diagnosis, treatment and prognostic prediction of HNSCC.

In the big data information era, the traditional matrix clustering tools, such as PCA and SVD, are not satisfactory and negative elements lack scientific explanation in application of practical problems. Clustering, which is divided into one-way or two-way clustering, is based on the principle that genes with comparable expression patterns have similar or related functions. It is one of the most important methods for processing gene expression data [5]. One-way clustering is whereby only rows or columns are clustered and its results are more influenced by unrelated columns or rows. Some of the commonly used one-way clustering algorithms include systematic clustering, self-organizing mapping clustering and principal component clustering. Two-way clustering is whereby the optimal set of sub-matrices are found in a matrix where rows and columns are significantly correlated. It allows overlap between classes, which is significant for gene chip data. Usually, a gene is not involved in only one biological process, it may be involved in multiple biological processes at the same time. Therefore, bidirectional clustering is more suitable for processing gene expression data. Non-negative matrix factorization is a two-way clustering process [6]. Compared to the other standard decomposition methods, non-negative matrix factorization (NMF) has 3 main advantages, namely, no parameters, good interpretability and good numerical results [6]. Based on gene expression profile data, non-negative matrix factorization has been widely used for cancer classification [7].

We performed molecular clustering and prognostic modeling of HNSCC samples from TCGA database and two validation groups (collected at the Department of Oral and Maxillofacial Head and Neck Oncology, Shanghai Ninth People’s Hospital and The First Affiliated Hospital of Zhengzhou University) based on NMF. This was aimed at appropriately informing the classification of HNSCC patients for treatment selection and prognostic prediction.

Data acquisition

RNA sequencing data, together with clinical and survival information of HNSCC patients were obtained from the TCGA Data Portal (https://portal.gdc.cancer.gov/repository). Post-operative tumor tissues and normal tissues were collected from 80 HNSCC patients from October 2009 to October 2016. Sixty patients diagnosed with HNSCC between 2015 and 2019 were collected from the First Affiliated Hospital of Zhengzhou University. Clinical information of the patients is shown in Tables 1, 2, and 3. Since this was a retrospective study, the informed consent requirement was waived.

Table 1

Basic clinical characteristics of derivation cohort.

Characteristic	Levels	Overall
N		502
T stage, n (%)	T1	33 (6.8%)
	T2	144 (29.6%)
	T3	131 (26.9%)
	T4	179 (36.8%)
N stage, n (%)	N0	239 (49.8%)
	N1	80 (16.7%)
	N2	154 (32.1%)
	N3	7 (1.5%)
M stage, n (%)	M0	472 (99%)
	M1	5 (1%)
Clinical stage, n (%)	Stage I	19 (3.9%)
	Stage II	95 (19.5%)
	Stage III	102 (20.9%)
	Stage IV	272 (55.7%)
Gender, n (%)	Female	134 (26.7%)
	Male	368 (73.3%)
Age, n (%)	≤60	245 (48.9%)
	>60	256 (51.1%)
Race, n (%)	Asian	10 (2.1%)
	Black or African American	47 (9.7%)
	White	428 (88.2%)
Age, median (IQR)		61 (53, 69)

Table 2

Basic clinical characteristics of validation cohort.

Characteristic	levels	Overall
N		80
Pathologic T stage, n (%)	T2	14 (17.5%)
	T3	32 (40%)
	T4	34 (42.5%)
Pathologic N stage, n (%)	N0	40 (50%)
	N1	16 (20%)
	N2	20 (25%)
	N3	4 (5%)
Pathologic M stage, n (%)	M0	51 (65.4%)
	M1	27 (34.6%)
Pathologic stage, n (%)	Stage II	39 (49.4%)
	Stage III	36 (45.6%)
	Stage IV	4 (5.1%)
Gender, n (%)	Female	35 (43.8%)
	Male	45 (56.2%)
Age, n (%)	≤60	40 (50%)
	>60	40 (50%)
Age, median (IQR)		61.5 (51, 74.25)

Table 3

Basic clinical characteristics of validation cohort from the First Affiliated Hospital of Zhengzhou University.

Characteristic	Levels	Overall
N		60
Pathologic T stage, n (%)	T2	12 (20%)
	T3	21 (35%)
	T4	27 (45%)
Pathologic N stage, n (%)	N0	27 (45%)
	N1	9 (15%)
	N2	12 (20%)
	N3	12 (20%)
Pathologic M stage, n (%)	M0	42 (70%)
	M1	18 (30%)
Pathologic stage, n (%)	Stage II	30 (50%)
	Stage III	21 (35%)
	Stage IV	9 (15%)
Gender, n (%)	Female	27 (45%)
	Male	33 (55%)
Age, n (%)	≤60	24 (40%)
	>60	36 (60%)
Age, median (IQR)		59.5 (50, 76.2)

Consensus clustering of HNSCC samples based on the NMF model

The NMF cluster was constructed using the Consensus Cluster Plus package [8]. Non-negative matrix factorization hierarchical clustering was performed using the adjusted and unified dataset, the number of clusters k values were from 2 to 8. Based on the clustering effect, the value with better clustering stability was selected [9]. With regards to NMF classification results, Kaplan-Meier survival analysis was performed. Differences in survival outcomes among different groups of patients with different immune cell infiltration levels were evaluated using the vioplot package in R.

Construction of the prognostic model

Differentially expressed genes (DEGs) were analyzed using edgeR, where they were screened using threshold values set to the absolute value of logFC >1 and FDR <0.05. DEGs that were significantly associated with overall survival (OS) outcomes in HNSCC patients were screened using univariate Cox regression analysis. Collinearity between genes was eliminated by Lasso regression analysis. Then, genes were included in a multifactorial Cox regression analysis model for further screening to identify predictive model component genes.

The prognostic signature was used as the risk score = .Where n is the number of prognostic genes, expi is the expression value of gene i, while βi is the regression coefficient of gene i in Cox regression analysis. The risk score was determined for every patient according to the formula. The median of the scores was the cut-off value, from which all patients were divided into high-risk and low-risk groups. Overall survival curves for the different groups of patients were plotted using the Kaplan-Meier method after which the log-rank test was performed.

The predictive ability of the proposed model was assessed using ROC curves and calibration plots. The previously described risk calculation formula was also used to calculate the risk score for every patient in the validation group. The ROC and calibration plots were also used to validate the predictive ability of the model.

HNSCC tissues were freshly isolated from surgical samples, and HNSCC diagnosis confirmed by pathology. Approximately 100 mg of samples from the tumor center were stored in liquid nitrogen, and paired with approximately 100 mg of normal tissue (>5 cm from the tumor tissue) samples from the same patient. Total RNA was extracted using the TRIzol method, and RNA concentrations in each sample measured by the Nano Drop 2000 system. Then, qRT-PCR was performed according to FastStart Universal SYBR Green Master operating instructions [10]. Beta-actin was the internal reference. Data were processed using the 2-ΔΔCT method.

Correlation analysis of model-independent prognostic and clinical characteristics

Univariate and multifactor Cox regression analyses of risk scores were performed to determine whether the model had an independent prognostic value. Incase the risk score was significantly different from OS in both univariate and multivariate Cox analyses, it was considered to be an independent risk factor. Finally, DCA was used to prove the clinical validity of the established model.

Immune cell infiltration analysis

Single-sample gene set enrichment analysis (ssGSEA) of a set of 16 immune-related genes was performed to quantify the activities and enrichment levels of immune cells, functions or related pathways in HNSCC. Expression analysis was performed to determine the association between the risk score and immune-related genes, such as m6a, ferroptosis, cellular autophagy, tumor mutation burden (TMB) and major histocompatibility complex (MHC). Based on the IMvigor210 immunotherapy cohort, which consisted of patients administered with the anti-PD-L1 antibody, Atezolizumab, we assessed the robustness of the classification and the ability to predict immunotherapeutic responses.

Results

Clustering based on the NMF model divided the samples into 2 subgroups

To reduce the impact of multicenter source and batch processing of samples, data were calibrated using “ComBat” in R [11]. With regards to clustering stability [12, 13], stability was found to be better when k = 2, therefore, k = 2 was used for judgment (Figure 1a). Survival curves and log-rank test results revealed that prognostic outcomes for the 2 subgroups were significantly different ((p < 0.05; Figure 1b-c). Moreover, immune cell infiltrations between the two subgroups were significantly different (Figure S1). Immune cell infiltration levels, including T, NK and CD8 cells, were higher in group C1, relative to C2.

Figure 1

(a) Non-negative matrix factorization cluster analysis. The best fitted cluster was k = 2 value. KM curves showing PFS (b) and OS (c) for 2 clusters. d. Tenfold cross-validated error (first vertical line equals the minimum error, whereas the second vertical line shows the cross-validated error within 1 standard error of the minimum) (left). The profile of coefficients in the model at varying levels of penalization plotted against the log (lambda) sequence (right).

Prognostic models

The above classification confirms differences in prognostic outcomes between the two clusters, therefore, DEGs between the two clusters were subjected to univariate Cox analysis to obtain prognosis-associated DEGs. LASSO was also used to screen the 13 associated genes (Figure 1d). These genes were subjected to multifactorial Cox analysis from which 9 DEGs and their correlation coefficients were obtained (Table 4). The prognostic model risk score was: risk score = 0.78∗ expression levels of HAUS6 -0.38∗ expression levels of SCNN1D+ 0.76∗ expression levels of S100A1 -1.64∗ expression level of TNFRSF4 -1.51∗ expression levels of FBX O 17 + 0.75∗ expression levels of IRF9+ 0.65∗ expression levels of IFI6 -0.41∗ expression levels of PTGS2+ 0.32∗ expression levels of MSC. The risk score for each patient was calculated based on the regression coefficients according to the prognostic model. Then, patients were assigned into high- and low-risk groups using median risk scores.

Figure 4

a. The relationships among tumor mutation burden, immune infiltration, and risk score. b. Immunotherapeutic responses of the high- and low-risk groups. c. Gene set enrichment analysis (GSEA, www.broadinstitute.org/gsea/). d. Functional network enrichment analysis.

Time-dependent ROC curves showed 1-, 3- and 5- year AUCs of 0.852, 0.890 and 0.953, respectively (Figure 2a). In this model, OS time of high-risk group patients was significantly shorter, compared to low-risk group patients (Figure 2b). A satisfactory agreement between the observed values was observed in the calibration curves (Figure 2c). Applying the same prognostic score to the validation set, Kaplan-Meier survival curves revealed that patients with high risk scores had lower OS, compared to those with low risk scores, and OS outcomes between the two groups were significantly different (Figure 3a). The 1- and 5-year AUC values for the validation set ranged from 0.767 to 0.862, indicating that the model had a good predictive performance in the external validation set (Figure 3b). Similarly, in the validation cohort from the Zheng University Hospital, the area under ROC curve for 3-year and 5-year survival rates were 0.766 and 0.765, respectively (Figure S2).

Figure 2

Prognostic analysis of the model in the derivation cohort. a. AUC of time-dependent ROC curves verified the prognostic performance of the risk score in the derivation cohort. b. Kaplan-Meier curves for OS of patients in the high-risk group and low-risk group in the derivation cohort. c. Calibration plot for model.

Figure 3

Validation of the model in the validation cohort. a. Kaplan-Meier curves for OS of patients in the high-risk and low-risk groups of the validation cohort. b. AUCs of time-dependent ROC curves verified the prognostic performance of the risk score in the validation cohort. c. Results of qRT-PCR analysis. d. The decision curve analyses (DCA) for clinical significance of this model.

PCR analysis showed that 9 genes were differentially expressed in the validation group (Figure 3c), in line with findings from the TCGA cohort. Findings from Kaplan-Meier survival curve analyses for the 9 genes are shown in Figure S3. These findings suggest that the risk score model has good sensitivity and specificity for prognostic prediction of HNSCC. Multifactorial Cox regression analysis was performed by combining clinical indicators of the patients (risk score, age, gender, stage and grade among others). The risk score was associated with survival outcomes (Table 5). Then, decision curves were used to determine the clinical net benefit of the model. Decision curve analysis showed that the model was clinically useful (Figure 3d). In conclusion, independent of other clinical factors, the risk score is a potential prognostic indicator for HNSCC.

Table 5

Univariate and multivariate Cox regression models were used to detect the prognostic elements.

Characteristics	Univariate analysis		Multivariate analysis
	Hazard ratio (95% CI)	P value	Hazard ratio (95% CI)	P value
T stage
T1	Reference
T2	1.086 (0.568–2.074)	0.803
T3	1.461 (0.769–2.773)	0.247
T4	1.249 (0.665–2.344)	0.490
N stage
N0	Reference
N1	1.058 (0.728–1.539)	0.768	0.999 (0.682–1.465)	0.997
N2&N3	1.404 (1.038–1.900)	0.028	1.469 (1.077–2.003)	0.015
M stage
M0	Reference
M1	4.745 (1.748–12.883)	0.002	4.288 (1.563–11.761)	0.005
Age
≤60	Reference
>60	1.252 (0.956–1.639)	0.102
Gender
Female	Reference
Male	0.764 (0.574–1.018)	0.066	0.779 (0.579–1.046)	0.097
Riskscore (low vs high)	0.770 (0.672–0.883)	<0.001	0.757 (0.660–0.870)	<0.001
Clinical stage
Stage I&Stage II	Reference
Stage III&Stage IV	1.217 (0.878–1.688)	0.238

Immunogenesis and enrichment analysis

Analysis of the relationship between the risk score and m6a, ferroptosis, cellular autophagy as well as other related genes revealed that the risk score was closely associated with immune-related genes (Figure S4). The TMB refers to the number of base mutations per million bases and is a marker for the efficacy of immune checkpoint inhibitors. The higher the TMB, the more neoantigens can be recognized by T cells and the better the immunotherapeutic effect. We found a negative correlation between the risk score and TMB, which may explain the poor prognostic outcomes for high-risk patients (Figure 4a). There was a higher probability of higher benefit for high-risk patients subjected to immunotherapy (complete response (CR), partial response (PR), no clinical benefit (progressive disease (PD) or Stable Disease (SD)). This provides new options for treatment of patients with subsequent tumors (Figure 4b). GSEA showed that the high-risk group was enriched in dilated cardiomyopathy whereas the low-risk group was mainly associated with tumorigenesis. These findings may partially explain the biological differences between the low- and high-risk groups at the genetic and pathway level (Figure 4c). Enrichment and signaling pathway analyses were performed for DEGs to elucidate on their biological significance. They were found to be mainly enriched in cell cycle, DNA replication, catalytic activities, acting on DNA, chromosomal region, human papillomavirus infection, organelle fission and condensed chromosome (Figure 4d). ssGSEA showed that the high-risk group had higher levels of infiltrating immune cells, especially T helper cells, macrophages, regulatory T (Tregs) cells and tumor-infiltrating lymphocytes (TILs) (Figure 5a). In the TCGA cohort, apart from MHC_class_I and type 1 IFN response pathways, activities of the other 11 immune pathways in the high-risk group were lower than those of the low-risk group (Figure 5b).

Figure 5

Comparisons of ssGSEA scores between different risk groups in the derivation cohort. Scores of 16 immune cells (a) and 13 immune-related functions (b) are shown in boxplots. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ns = not significant.

Drug sensitivity analysis

The highest negative correlation score was for chrysin (−0.776). Chrysin is a drug with various pharmacological activities, including anti-tumor, anti-inflammatory, anti-bacterial, anti-anxiety and anti-oxidant effects [14], suggesting a possible therapeutic effect in HNSCC. The next highest score was MS-275. Previous studies have shown that MS-275 has a selective killing effect on gastric adenocarcinoma cells [15], 1, 4-chrysenequinone (an Ahr-activator) and piperlongumine (inhibits tumor autophagy leading to reduced cell proliferation viability) (Table S1).

Discussion

Globally, HNSCC is a common malignancy, with 550,000 new cases and about 380,000 deaths per year [16, 17]. It is aggressive, lethal and causes serious facial deformities, speech, chewing and swallowing dysfunctions as well as psychosocial problems. Although radical surgical techniques, repair and reconstruction techniques for HNSCC have become increasingly sophisticated in the last 20 years, there have been no significant improvements in 5-year survival outcomes [18].

The prediction of individual patient prognosis will greatly inform treatment decisions. Based on the NMF model, we staged HNSCC patients into two subgroups. There were significant differences in OS outcomes between the two subgroups, with patients in the subgroup with more abundant immune cell infiltrations exhibiting better prognostic outcomes. A prognostic risk model consisting of nine genes was constructed. Patient scores were calculated based on the risk model and divided into two groups: high and low risk groups. There were significant differences in prognostic outcomes of patients in the two groups, with the prognostic outcomes of high-risk patients being significantly lower than those of low-risk patients. Moreover, ROC and calibration curves of the model achieved remarkable results, which revealed that the model has better discriminatory abilities. The DCA also showed that reliability and accuracy of the prediction model was better than that of the other clinical indicators.

TNFRSF4, one of the component genes of the model, is predominantly inducibly expressed in activated CD4+ and CD8+ cells [19]. Binding of TNFRSF4 to ligands promotes clonal proliferation of T cells, enhances T cell memory, proliferation, immune surveillance and killer cell expansion. However, it inhibits immune tolerance development [20]. In addition, TNFRSF4 expressing positive T cells can reduce suppressive factors in the tumor immune microenvironment and effectively inhibit tumor invasion and metastasis [21]. Expressions of TNFRSF4 in breast cancer, melanoma, and lymphoma have been discussed in previous studies [22, 23]. Targeting TNFRSF4 has a role in anti-breast cancer and melanoma treatment [24]. In glioma, hepatocellular carcinoma and lung adenocarcinoma, FBXO17 promotes cell proliferation, migration and invasion through the Akt/GSK-3β/Snail pathway [25, 26, 27]. Not only is IRF9 important for antiviral responses, it is also involved in autoimmunity [28]. IFI6, which promotes the metastatic potential of breast cancer cells through mtROS, belongs to the ISG12 gene family, which is composed of four members, ISG12a, ISG12b, ISG12c and IFI6 [29].

Immune cell infiltrations in tumor sites is the basis for effective immunotherapy [30]. Therefore, understanding immune cell infiltrations in the TME is key to improving response rates and developing new immunotherapeutic strategies [31]. Although T cell properties have been widely evaluated, other immune cells of the innate and adaptive immune system, including dendritic cells, macrophages, natural killer cells, and B cells also influence tumor progression and immunotherapeutic responses [31]. Elevated macrophage levels are associated with poor cancer prognosis [32]. Macrophage infiltrations in the tumor microenvironment promote tumor growth, angiogenesis, invasion and metastasis [33]. Due to their potent tumor-killing abilities, T cells are the focus of tumor immunity. Within the tumor microenvironment, different types of T cells, including cytotoxic T cells (CTL), T follicular helper cells (Tfh) and regulatory T cells (Tregs) are involved in T cell-mediated immune responses [34]. Tumor-infiltrating lymphocytes are positively associated with survival outcomes in various cancers; however, due to immunosuppression of the tumor microenvironment, tumor-infiltrating T cells are often unable to control tumor growth, leading to their depletion or dysfunction [35, 36, 37]. Enrichment analysis revealed that the DEGs are mainly enriched in the cell cycle, DNA replication, catalytic activities and acting on DNA while GSEA showed that low-risk patients were predominantly enriched with immunodeficiency and tumor-associated pathways.

The risk scores were strongly associated with m6a-related, ferroptosis-related and autophagy-related genes. Studies on tumor immunotherapy, which is an important area of research, have rapidly progressed [38]. Immune-suppressants such as PD-1/PD-L1 have successfully been developed [39]. In this study, responsiveness to PD-1/PD-L1 in both high and low risk patient groups revealed that high risk patients responded to immunotherapy better than low risk patients.

Conclusions

Based on the NMF algorithm, we screened for DEGs and constructed an associated prognostic model that can independently predict prognostic outcomes for HNSCC patients. The predictive performance of the model was found to be stable and could inform individualized treatment of HNSCC patients. Furthermore, the genes in the prognostic risk model are potential immunotherapeutic targets for HNSCC.

Declarations

Author contribution statement

Xin-Yu Li; Hong-Bang An; Xi-tao Yang: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Lu-yu Zhang; Hui Liu; Yu-chen Shen: Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement

Data will be made available on request.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.

Table 4

Gene correspondence coefficient.

id	coef
HAUS6	0.781261442
SCNN1D	-0.382482986
S100A1	0.760949529
TNFRSF4	-1.642657948
FBXO17	-1.512965493
IRF9	0.751127416
IFI6	0.654360943
PTGS2	-0.41098144
MSC	0.322998811