Overexpressed PLAU and its potential prognostic value in head and neck squamous cell carcinoma

Data collection and normalization

GSE25099 from 79 HNSCC patients (Peng et al., 2011), GSE13601 consisting of 37 HNSCC patients (Estilo et al., 2009), GSE65858 from 270 HNSCC patients (Wichmann et al., 2015), GSE136037 from 49 HNSCC patients (Alfieri et al., 2020), and The Cancer Genome Atlas (TCGA) -HNSC cohorts including 546 HNSCC HTSeq-counts, methylation profiles and related clinical information were downloaded. CalcNormFactors was used to calculate normalization factors to scale the gene expression in TCGA dataset (Anders & Huber, 2010; Le et al., 2019). Youden index (sensitivity + specificity -1) was used to calculate the best cutoff of survival analysis by R package (“SurvivalROC”) (Huang, Liao & Li, 2017; Luo & Xiong, 2013). The GSE65858 and GSE136037 datasets from the GEO database were converted to transcripts per million (TPM) (Zhao, Ye & Stanton, 2020a). And the term “N-” means primary lesions in HNSCC patients without neck lymph node metastasis, and “N+” means primary lesions in HNSCC patients with neck lymph node metastasis. The normalized data provided from original studies in other datasets were used in this study directly.

Survival analysis

Briefly, the expression of PLAU was categorized into low or high by ‘Maximally Selected Rank Statistics’ (maxstat) method (Lausen & Schumacher, 1992). The cutoff value of PLAU mRNA expression in TCGA was 9,532 in Table 1. Older or younger is classified based on its mean age. The overall survival (OS), progression-free interval (PFI) or recurrence free survival (RSF) curves were visualized by Kaplan–Meier plots. Univariate and multivariate Cox regression analysis was applied to death hazard ratios calculation after the proportional hazard assumption was tested.

Go analysis and co-expression network establishment

Go enrichment analysis was performed by Bioconductor package “clusterProfiler” (Yu et al., 2012). The method of KEGG enrichment analysis was performed same as the Go enrichment analysis. The co-expression genes were screened using R packages (“limma”).From the TCGA-HNSCC databases, we used Pearson correlation coefficients (—Pearson correlation coefficient— > 0.5 and P-value <0.001) and the z-test to examine the correlation between PLAU expression level and co-expression genes. STRING database and Cytoscape tool were used to construct the protein-protein-interactions (PPI) of genes (Struk & Jacobs, 2019).

Immune cell environment analysis

Estimation of stromal and immune cells in malignant tumor tissues using Expression data (ESTIMATE) is a tool for predicting tumor purity, and the presence of infiltrating stromal/immune cells in tumor tissues using gene expression data (Yoshihara et al., 2013). We used R packages (“estimate”) to get immune-environment scores in HNSCC patients. R packages (“CIBERSORT.R”, https://cibersortx.stanford.edu/) (Newman et al., 2015) was used to deconvoluted 22 common immune cell proportions in HNSCC patients. The correlation of PLAU expression with immune environment and 22 immune cells were investigated by Pearson test.

Software and statistical analyses

GraphPad Prism 8 or R studio (version 3.5.3) was used to evaluate all data (Kamdem et al., 2019; Le & Huynh, 2019). Chi-squre or Fisher’s exact test was performed to compare the differences expression of PLAU across different groups. P < 0.05 was considered statistically significant. The detailed codes and other packages’ information were included in the supplementary materials.

PLAU mRNA is over-expressed in HNSCC

To discover the expression of PLAU mRNA in HNSCC, we analyzed three independent patient cohorts, which showed a consistent result that PLAU mRNA expression was elevated in HNSCC tumors than normal tissues(all p < 0.01, Figs. 1A–1C). And the overexpression of PLAU mRNA was confirmed in 10 different HNSCC cell lines as compared to 4 types of human keratinocyte cell lines (p = 0.002, Fig. 1D). Furthermore, PLAU mRNA in HPV positive HNSCC samples was interposed between adjacent normal tissues and HPV negative tumors (Fig. S1A).

There was significantly less PLAU mRNA in HPV positive tumors than HPV negative ones according to another patient cohort (Fig. S1B).

Association of PLAU mRNA with neck lymph node status in HNSCC

In order to study the role of PLAU mRNA may play in HNSCC patients, the relationship of PLAU mRNA and clinical parameters was further characterized in TCGA-HNSC cohort. As shown in Table 2, there was no significantly different expression of PLAU mRNA in different age, gender, clinical stage, and tumor stage, but the difference of PLAU in HPV positivity and neck node status was considered as significant (p = 0.001, p = 0.033, respectively). Higher expression of PLAU was founded in patients with neck lymph node metastasis than patient without neck lymph node metastasis in another independent cohort as well (Fig. 1E).

PLAU is an independent predictor of HNSCC prognosis

Considering the above findings, we continued to analysis the possible correlation of HNSCC survival time and PLAU expression. As shown in Table 3, age, clinical stage, tumor size, neck lymph node status, HPV positivity and PLAU expression was considered to be significantly associated with overall survival time in univariate cox analysis of 496 HNSCC patients. Multivariate cox analysis indicated the hazard ratio of death was reached to 1.52 when high PLAU expressed HNSCC patients compared to patients with low PLAU expression after excluding the potential affections from age, tumor size, neck node metastasis and HPV positivity (Fig. 2A, p = 0.012, 95% CI [1.09–2.10]).

The survival curves of PLAU expression were visualized by Kaplan–Meier plots (Fig. 2B), which implied that high PLAU expressed HNSCC patients had decreased overall survival probability than patients with low PLAU expression (Fig. 2B). The same finding could be observed in another independent HNSCC cohort (Fig. 2C). Additionally, high expression of PLAU in HNSCC patients was founded to predict unfavorable outcomes in terms of PFI and RSF (Figs. 2D–2E).

Performance of PLAU expression in predicting 5-year overall survival outcomes of HNSCC patients

The predictive performance of PLAU expression in 5-year overall survival outcomes of HNSCC patients was analyzed by ROC analysis. As shown in Fig. 3, the areas under the ROC curve (AUCs) of PLAU expression was 0.795, higher than HPV status (AUC: 0.531), neck node status (AUC: 0.58) and tumor size stage (AUC: 0.532). Combining PLAU expression and neck node status as an integrated factor (N+PLAU), the AUC was 0.862, with an 88.5% sensitivity, 64% specificity, 70.4% positive predictive value and 80.1% negative predictive value. Taken together, PLAU itself or combined with neck node status may facilitate the overall survival prediction for HNSCC patients.

Hypomethylation of PLAU in HNSCC patients

A growing body of literature reports that DNA methylation and mRNA expression most likely to be negative correlated (Huang et al., 2015). Methylation profiles of 496 TCGA-HNSC patients were analyzed to investigate the possible upstream reasons for elevated PLAU mRNA expression in HNSCC tumors. As shown in Fig. 4, the methylated level of PLAU was successively and significantly decreased in adjacent normal tissues, HPV positive tumors and HPV negative tumors (Fig. 4A). The PLAU methylation and mRNA expression was negatively correlated (Fig. 4B, R =  − 0.43, p < 0.001). As expected, patients with hypomethylated PLAU, indicating higher expression of PLAU mRNA, had worse overall survival outcomes than patients with hypermethylated PLAU (Fig. 4C).

Network establishment for PLAU correlated genes in HNSCC

To further understand the possible downstream reasons a total of 205 genes was correlated with PLAU expression in TCGA-HNSCC patients (21 of 205 genes were negatively correlated with PLAU and 184 genes were positively correlated with PLAU). The top 20 genes of positively or negatively correlated with PLAU are shown in a Heatmap (Fig. 5A). The interaction network of these 205 genes was established based on STRING and Cytoscape (Fig. 5B). Next, we performed GO and KEGG enrichment analysis to understand the potential biological functions of PLAU in HNSCC. GO analysis showed that the major biological process (extracellular matrix organization), cellular component (extracellular matrix), and molecular functions (cell adhesion molecule binding) may contribute to PLAU related biology (Fig. 5C). KEGG pathway analysis illuminated that the P13K-Akt signaling pathway, human papillomavirus infection, proteoglycans in cancer, and focal adhesion as significantly enriched by the PLAU co-expressed genes (Fig. 5D).

Distributions of tumor infiltrating immune cell in HNSCC patients with different PLAU expression

More and more evidence revealed the tumor immune microenvironment is a crucial factor in tumor biology (Mao et al., 2020). To interpret the role of PLAU expression in HNSCC based on immunity conception, the scores or proportions of tumor infiltrating cells were compared in TCGA-HNSCC cohort. It shows that PLAU expression was positively correlated with stromal score (Fig. 6A). Further analysis found the expression level of PLAU was positively correlated with M1 type macrophages, negatively associated with CD4 + T cell, Tregs cell, and follicular helper T cell (Figs. 6B–6E) (All p < 0.05).