A novel focal adhesion related gene signature for prognostic prediction in hepatocellular carcinoma

Hepatocellular carcinoma (HCC) is a highly heterogeneous disease. Reduced expression of focal adhesion is considered as an important prerequisite for tumor cell invasion and metastasis. However, the prognostic value of focal adhesion related genes in HCC remains to be further determined. In this study, RNA expression profiles were downloaded from public databases. A five focal adhesion related gene signature model was established by the least absolute shrinkage and selection operator Cox regression analysis, which categorized patients into high- and low-risk groups. Multivariate Cox regression analysis showed that the risk score was an independent predictor for overall survival. Single-sample gene set enrichment analysis revealed that immune status was different between the two risk groups, and tumor-related pathways were enriched in high-risk group. The risk score was significantly associated with tumor grade, tumor stage, immune scores, and immune infiltrate types. Pearson correlation showed that the expression level of prognostic genes was associated with anti-tumor drug sensitivity. Besides, the mRNA and protein expression of prognostic genes was significantly different between HCC tissues and adjacent non-tumorous tissues in our separate cohort. Taken together, a novel focal adhesion related gene signature can be used for prognostic prediction in HCC, which may be a therapeutic alternative.


INTRODUCTION
Liver cancer is the third global leading cause of cancerrelated death globally with a high mortality rate [1]. As the main histologic type of liver cancer, hepatocellular carcinoma (HCC) accounts for approximately 90% of all liver cancers and causing approximately 750,000 deaths annually [2,3]. In most developing countries, hepatitis B virus (HBV) and hepatitis C virus (HCV) infections are the main risk factors of HCC, while in most developed countries, alcohol abuse and metabolic syndrome-related diabetes and obesity are the main underlying disorders [4]. Surgical resection, liver transplantation, and radiofrequency ablation remain the mainstay of treatment for HCC [5,6]. Although an increasing number of diagnostic and therapeutic strategies have been developed, the prognosis of HCC patients remains poor due to the elusive molecular mechanisms. It is therefore necessary to gain further insights into the molecular mechanisms underlying HCC progression and seek novel prognostic biomarkers and therapeutic strategies for this fatal disease.
Focal adhesion, as the main connection between cells and the extracellular matrix (ECM), helps maintain the tension of cells during movement and signal transmission of cell survival [7,8]. Of various microenvironmental factors affecting cancer cell AGING resistance, cell adhesion to the ECM has recently been identified as a key determinant [9]. Reduced expression of focal adhesion promotes the process of epithelialmesenchymal transition and is considered as an important prerequisite for tumor cell invasion and metastasis [10][11][12][13][14]. Focal adhesion molecules including integrins, focal adhesion kinases, and growth factor receptors have been shown to be associated with cancer progression [15][16][17]. Besides, more and more research has demonstrated that focal adhesion related genes such as ITGA [18,19], NEK2 [20][21][22], and cPLA2α [23,24] all play a promoting role in HCC. However, whether the focal adhesion related genes are associated with the prognosis of HCC patients remains unclear. Given the most upstream localization, receptors within focal adhesions were supposed to be an ideal starting point for therapeutic intervention [9], suggesting that focal adhesion may prove to be the prognostic indicator and a novel therapeutic target in HCC patients.
In this study, we developed a prognostic signature for HCC patients based on focal adhesion related genes. We also investigated the correlation of the prognostic model risk score with the clinicopathological characteristics. In addition, we analyzed the role of prognostic model risk scores in immune infiltrate types, tumor stemness features, and prognostic gene expression in drug chemoresistance. Also, the stability and reliability of the model were demonstrated in an independent and external validation cohort, and the mRNA and protein expressions of the prognostic genes between HCC and adjacent non-tumorous tissues were verified by our laboratory experiments. It is our hope that the findings of the present study could provide a novel prognostic model and therapeutic targets for HCC.

RESULTS
The flow chart of this study is demonstrated in Figure 1. Altogether 365 HCC patients from The Cancer Genome Atlas hepatocellular carcinoma cohort (TCGA-LIHC) and 231 HCC patients from the International Cancer Genome Consortium hepatocellular carcinoma cohort (ICGC-LIRI-JP) were eventually recruited for analysis. The clinicopathological characteristics of these patients are detailed in Table 1.

Identification of DEGs in the TCGA cohort
More than 50% of focal adhesion related genes (131/199, 65.8%) were identified as differentially expressed genes (DEGs) between the tumor and adjacent non-tumorous tissues using the "limma" R package, and Univariate Cox regression analysis showed that 38 of the 131 DEGs were correlated with overall survival (OS) ( Figure 2A). As the expression level of ACTN3 was 0 in more than 80 samples, it was excluded from this analysis. Finally, 37 focal adhesion related DEGs were preserved ( Figure 2B-2C). The protein-protein interaction network and correlation network between these DEGs are shown in Figure 2D-2E.

Construction of a prognosis model in the TCGA cohort
The least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed to construct a prognosis model based on the expression profiles of the above-mentioned 37 DEGs. Using the optimal value of λ, a 5-gene signature was identified (Supplementary Figure 1). The risk score was calculated by using the following equation: Score = ˗0.031 * expression level of FYN + 0.170 * expression level of PPP1CB + 0.003 * expression level of PPP1CC + 0.084 * expression level of RAC1 + 0.045 * expression level of SPP1. The patients were categorized into a high-risk group and a low-risk group based on the median cut-off value ( Figure 3A). The tumor grades were higher and the TNM stages were more advanced in the TCGA cohort of high-risk group as compared with those in low-risk group (Table 2). PCA and t-SNE analyses showed that patients in these different risk groups were distributed in two directions ( Figure 3E-3F). As shown in Figure  3B, death usually occurred earlier in high-risk group than that in low-risk group (P = 0.007). Conformably, the Kaplan-Meier curve indicated that the prognosis of high-risk patients was significantly poorer than that of low-risk patients (P < 0.001) ( Figure 3I). The predictive performance of the risk score for OS was evaluated by time-dependent ROC curves. The area under the curve (AUC) of the prognosis model at 1, 2 and 3 years was 0.738, 0.681 and 0.652, respectively, indicating a good predictive performance ( Figure 3J). Survival analysis using the optimal cut-off expression value of each prognostic gene showed that high expressions of PPP1CB, PPP1CC, RAC1 and SPP1 were correlated with poor prognoses, whereas FYN showed an opposite correlation (P < 0.05) (Supplementary Figure 2A-2E). The different expression of each prognostic gene between HCC and adjacent non-tumorous tissues in the TCGA cohort is shown in Supplementary Figure 3.

Confirmation of the 5-gene signature in the ICGC cohort
The predictive stability of this prognosis model was further verified by using HCC samples from the ICGC database. Consistent with the above results, the HCC patients categorized into a high-risk group and a AGING low-risk group based on the median value from the TCGA cohort ( Figure 3C). The TNM stage was more advanced in in the ICGC cohort of high-risk group (Table 2). PCA and t-SNE analyses demonstrated that patients in the two subgroups were distributed discretely ( Figure 3G-3H). Patients in high-risk group died earlier ( Figure 3D) and had shorter survival durations ( Figure  3K) than patients in low-risk group (both P < 0.05). Additionally, the AUC of the 5-gene signature at 1, 2 and 3 years was 0.681, 0.661 and 0.673, respectively ( Figure 3L). The result of survival analysis of each prognostic gene with high or low expression is shown in Supplementary Figure 2F-2J.

Relationship between the clinical characteristics and risk score
The correlation between the risk score and various clinical characteristics in HCC patients based on TCGA and ICGC data is indicated in Figure 5A-5G, respectively. Higher risk score was observed in patients with the higher tumor grades (P < 0.001) and more advanced tumor stage (P = 0.004) in the TCGA cohort. However, the risk score had no significant correlation with age and gender. The same analysis in the ICGC cohort showed the similar results (There were no data about the grade of LICH in the ICGC cohort). Further analysis of the correlation of the expression level of the prognostic genes with age, gender, tumor stage and tumor grade of HCC patients showed that the expression of FYN, PPP1CC, RAC1, and SPP1 was significantly correlated with the tumor grade (P < 0.05). In addition, the expression of PPP1BC, PPP1CC and RAC1 was correlated with the tumor stage (P < 0.05) (Supplementary Figure 4C-4D); the expression of PPP1CC and SPP1 was correlated with age (P < 0.05) (Supplementary Figure 4A); the expression of PPP1CB was correlated with gender (P < 0.05) (Supplementary Figure 4B).

Analysis of the immune status and tumor microenvironment
To further explore the immune status in different risk groups, the enrichment scores of diverse immune cell subpopulations and related functions were quantitated by single-sample gene set enrichment analysis (ssGSEA). The ratios of aDCs, DCs, iDCs, macrophages, Th2 cells, and Treg were high in patients of high-risk group, and vise versa for the score of Mast cells (P < 0.05) ( Figure 6A-6B). Similarly, the score of immune-related functions such as APC coinhibition, APC co-stimulation, Check point, HLA, MHC class I, Parainflammation, T cell co-inhibition, and T cell co-stimulation was significantly higher in high-risk group, and vise versa for the score of Type II IFN (P < 0.05) ( Figure 6C-6D). All these results suggest that the immune status and tumor microenvironment might contribute to the prognosis of HCC patients with high expressions of focal adhesion related genes.
Subsequently, we examined the risk score distribution in different immune subtypes in HCC patients. Of the six types of immune infiltrates identified in human tumors [C1 (wound healing), C2 (INF-γ dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-β dominant)] [25], and no patient sample belonged to the C5 immune subtype in HCC. The value of the risk score from C1 to C6 subtype decreased gradually. As shown in Figure   6E, the high-risk score was significantly correlated with C1 and C2, and the low-risk score was correlated with C3, C4 and C6 (P < 0.05). Furthermore, the high levels of PPP1CB, PPP1CC, RAC1, and SPP1 were positively correlated with C1, suggesting their promoting roles in HCC. On the contrary, the high level of FYN was negatively correlated with the C1 (P < 0.05) (Supplementary Figure 5), suggesting its tumorsuppressive role in HCC.
To explore whether the risk score was related to tumor stem cells and the immune microenvironment, the result of correlation analysis was showed that the risk score was positively correlated with immune scores, suggesting an association of the prognostic model with immune composition (P < 0.001). On the contrary, the risk score had no significant correlation with stemness score based on mRNA expression (RNAss), DNA methylation based stemness score (DNAss), and stromal score ( Figure 6F). Moreover, the correlations of tumor stem cells and the immune microenvironment with prognostic gene expression were analyzed. It was found that FYN, PPP1CB, and PPP1CC were correlated with RNAss, while FYN, PPP1CC, and SPP1 were correlated with the stromal score (P < 0.05), where FYN exhibited strongest correlation (R = 0.55) (Supplementary Figure 6). Furthermore, FYN, RAC1, and SPP1 were found to be positively correlated with the immune score which measures the presence of infiltrating immune cells (P < 0.05) (Supplementary Figure 6).

Analysis of biological functions and pathways
Using gene set enrichment analysis (GSEA), gene ontology (GO) enrichment and Kyoto encyclopedia of genes and genome (KEGG) pathways were analyzed in high-and low-risk groups. It was found that 20 main GO functions were enriched in high-risk group with a false discovery rate <0.05. The enriched GO functions showed that the risk score was significantly associated with the activations of regulation of autophagic regulation, cell cycle phase transition, intrinsic apoptotic signaling pathway, and positive regulation of Wnt signaling pathway ( Figure 7A). As shown by KEGG pathway analysis, the pathway enrichment was correlated with the cancer procession including cell cycle, MAPK signaling pathway, mismatch repair, mTOR signaling pathway, notch signaling pathway, P53 signaling pathway, pathways in cancer, VEGF signaling pathway, and Wnt signaling pathway. Of note, the focal adhesion pathway was also enriched in highrisk group ( Figure 7B). These results indicate that the prognostic model had an extensive influence of on the global transcriptome of HCC tissues.

Chemotherapy sensitivity analysis
Then we investigated the association between focal adhesion related genes and chemotherapy sensitivity in the NCI-60 database. The results showed that all the prognostic genes were to some extent correlated with chemotherapy drug sensitivity (Figure 8, Supplementary

Verification of the expression of prognostic genes between HCC and adjacent non-tumorous tissues
To validate the different expression of the five prognostic genes (FYN, PPP1CB, PPP1CC, RAC1, and SPP1) between HCC and adjacent non-tumorous tissues, real-time quantitative-polymerase chain reaction (qRT-PCR) and immunohistochemistry (IHC) were performed to analyze the mRNA and protein expression, respectively. As shown in Figure 9A, PPP1CB, PPP1CC, RAC1, and SPP1 were highly expressed in HCC tissues vs. adjacent non-tumorous tissues, while FYN was at a low level in HCC tissues. Immunohistochemical staining showed the same results ( Figure 9B). The validation results were consistent with the RNA sequencing results of the five prognostic genes in the TCGA cohort (Supplementary Figure 3).

DISCUSSION
HCC is one of the most frequent cancer type worldwide [26]. Although multiple molecular mechanisms have been found to be involved in HCC progression [27], the underlying mechanism remains unclear. Our study provided a systematical analysis of the different expression of 199 focal adhesion related genes between HCC and adjacent non-tumorous tissues and their relationships with OS. Then, a 5-gene signature model was constructed by the LASSO algorithm. Both the internal and external validation cohorts demonstrated that the model worked stably and exhibited consistent predictive performance. More importantly, it could be used as an independent risk factor for predicting the prognosis of HCC patients.
Among the prognostic genes, PPP1CB, PPP1CC, RAC1, and SPP1 were found as risk factors, and FYN was found as a protective factor. FYN is a member of the Src family kinases, which is involved in the focal adhesion signaling pathway. It was reported to be related to numerous solid tumors, and increased expression of FYN was found to play a promoting role both in cancer occurrence and progression, as well as be a mechanism of resistance to anticancer agents [28] such as breast cancer [29], and pancreatic cancer [30]. PPP1CB is a regulator of endothelial cell migration and plays a critical role in the angiogenic process. Lacobazzi et al. reported that PPP1CB inhibition could inhibit endothelial cell migration through focal adhesion turnover via the actin polymerization pathway [31]. Besides, recent studies have demonstrated that PPP1CC encodes protein phosphatase 1c, which may promote cancer cell proliferation through activating mutations in p53 [32]. RAC1 is identified as a master regulator of cell migration and anchorage-independent growth [33]. Earlier studies reported that SPP1 could promote cancer progression through modulating the expression of VEGF and regulating ECM protein [34,35]. Interestingly, all the five genes that we analyzed above were all enriched in the focal adhesion signaling pathway which is the molecular bridge mediating two-way crosstalk between ECM and cytoskeleton [36]. But other studies argued that focal adhesion signaling pathway promoted tumor progression and metastasis [37,38]. Therefore, we subjected the 5-gene signature model to predict the prognosis for HCC, anticipating that it could become a target for cancer drug therapy. We demonstrated that the 5-gene signature model could add prognostic value to the clinicopathological prognostic characteristics. GSEA enrichment analysis suggested the existence of potential biological processes and pathways in high-risk group. Consistently, our study identified the enrichment in the high-risk-score group including the MAPK, notch, P53, VEGF and Wnt signaling pathways. All these pathways are classical signaling pathways that have been identified to be involved in cancer process [39][40][41]. The Wnt signaling pathway is severely dysregulated in solid and non-solid malignant tumors, acting as part of the proliferation pathway that characterizes cancer cells [42]. Pérez-Plasencia et al. reported that the Wnt signaling pathway, as a common pathway in many cancer types, was a main attractive target in cancer therapy [43], including HCC [35]. Dysregulation of the MAPK cascade is linked to some key signaling components and phosphorylation events which can activate the process of tumorigenesis [44]. At least three different MAPK signal transduction pathways participate in the modulation and transduction of extracellular signals into the nucleus to induce response genes in such mammalian cells as ERK1/2, JNK1/2/3, and p38 [44,45], and all of them have proved to be correlated with HCC [46][47][48]. The Notch signaling pathway plays a key regulatory role in cellular fate, survival, and injury response in HCC cells [49,50]. Increasing evidence demonstrates that inhibition of the Notch signaling can enhance the therapeutic efficacy of antitumor agents in HCC cells [51,52]. As one of the most frequently mutated genes in HCC and other cancers, the tumor suppressor p53 modulates cellular stress conditions and responses to DNA damage and other cytotoxic stresses [53][54][55][56][57]. VEGF, one of the most prominent regulators involved in vasogenesis, was found to be highly expressed in human HCC specimens [58,59]. Interestingly, we found that the focal adhesion pathway was also enriched in HCC. Activation of the focal adhesion pathway mediated cancer cell survival, invasion, proliferation, and drug resistance [60][61][62].
And it has been shown to activate different cytoplasmic signaling pathways and co-regulate the pro-survival mechanisms [61,63,64]. Integrins are key mediators of cell adhesion molecules which control critical cell functions such as survival and migration by activating certain signaling mechanisms [65,66]. Given these findings, it is plausible to assume that focal adhesion has a close connection with HCC procession. However, more research is needed to clarify the specific role of the Focal Adhesion Pathway in HCC procession. Indeed, it is essential to further explore the efficacy of the combination of focal adhesion related genes with the long-term survival of HCC patients.
It is an urgent task to seek new biomarkers for early detection, prognostic assessment and decision making in the treatment of HCC. It was found that the abnormal expression of focal adhesion related prognostic genes was associated with poor prognosis of HCC patients [67], suggesting that these abnormally expressed genes may prove to be prognostic biomarkers in HCC. The tumor microenvironment and dysregulation of immune status also impact tumor progression. Increased numbers of studies suggest that the integration of tumor-infiltrating immune cells and clinicopathological characteristics may prove to be a potential prognosis model for predicting the response of patients to immune therapy [68]. In the present study, we found a positive correlation between the risk score and tumor-infiltrating immune cells suggesting that the prognostic genes for predicting the clinical outcomes of HCC patients may be associated with immune cell infiltration. Moreover, our data have also documented that the risk score is associated with the enrichment scores of diverse immune cell subpopulations and related functions, such as aDCs, DCs, iDCs, APC co-stimulation, HLA, and MHC class I, suggesting that the risk score may potentially affect the antigen presentation. Furthermore, we observed that the risk score was relevant to T cell subsets, like Th2 cells and Treg cells, indicating its potential role in regulating the secretion of cytokines from the Th cells.
Knowing that cancer stem cells (CSCs) are highly resistant to conventional chemotherapeutic drugs and radiation therapy, and promote cancer progression owing to their potent self-renewal and invasion capabilities [69,70], selection of CSCs with specific markers or signaling pathways could be an effective therapeutic strategy for the treatment of chemotherapyresistant HCC [71]. Of various microenvironmental factors affecting chemotherapy resistance, cell-cell adhesion and communication have recently been identified as key determinants [9,72]. In the present study, we investigated the expression of prognostic genes with stem-cell-like features measured by RNAss and DNAss. We found that FYN and PPP1CB were negatively correlated with RNAss, while PPP1CC was positively correlated with RNAss, which indicated that prognostic genes might have a relationship with cancer cell sensitivity or resistance to chemotherapy treatment. Moreover, our data showed that increased expression of the prognostic genes was associated with increased drug resistance or sensitivity to a number of chemotherapy drugs. Based on these observations, we preliminarily concluded that overexpression of the prognostic genes could decrease cancer cell sensitivity or resistance to drug treatment and be used as therapeutic targets.

CONCLUSIONS
In summary, the focal adhesion related gene signature identified in this study showed good performance in different cohorts and could be used to improve the AGING prediction accuracy of OS in HCC patients, though the action mechanism between the focal adhesion related genes and tumor immunity in HCC needs to be further explored. Taken together, would help gain deeper insights into the prognostic value and biological function of focal adhesion related genes in HCC and provide new possibilities for HCC therapeutic intervention.

Data collection
The RNA sequencing expression profile and corresponding clinical information were retrieved from two public cohorts: TCGA-LIHC and ICGC-LIRI-JP. A total of 365 HCC patients with liver cancer datasets were collected from the TCGA portal (https://portal.gdc.cancer.gov/repository).

Establishment of a focal adhesion related gene signature
Using the "limma" R package, DEGs between tumor and adjacent non-tumorous tissues were identified with a false discovery rate of <0.05 in the TCGA cohort. The prognosis genes in HCC patients were screened buy univariate Cox analysis. A protein-protein interaction network of the overlapping prognostic DEGs was constructed with similar interaction patterns grouped. The possibility of interactions between genes was further screened by LASSO algorithm. A prognostic model was constructed with the "glmnet" R package. The basic aim of LASSO was variable selection and some regression coefficients could be strictly equal to 0, thereby obtaining an interpretable model. The response variables were OS and the status of patients in the TCGA cohort, and the independent variables in the regression were the normalized expression matrix of candidate prognostic DEGs. The penalty parameter (λ) for the model following the minimum criteria (i.e. the value of λ corresponding to the lowest partial likelihood deviance) was determined by 10-fold cross-validation, a standard method to estimate the adjustment parameter λ. The risk score of the patients was calculated based on the expression level of each gene and its corresponding regression coefficient by using the following equation: Score = e sum (each gene's expression × corresponding coefficient) . The median risk score was used to classify patients into a high-risk group and a low-risk group. OS of the HCC patients was determined by survival analysis using the "survminer" R package. The predictive power of the model was assessed by time-dependent ROC curve analysis using the "survivalROC" R package. Based on the expression of the prognostic genes, we performed PCA and t-SNE with the "Rtsne" R package to explore the distribution of different groups. Furthermore, univariate and multivariate Cox regression analyses were performed to explore the independent prognostic value of the 5-gene signature.

Enrichment analysis
To further understand the immune function, ssGSEA was performed to calculate the score of immune cell infiltration and the activity of immune-related pathways between high-risk and low-risk groups using the "GSVA" R package. GO and KEGG analyses in highrisk and low-risk groups were performed using GSEA by GSEA software 4.1.

Tumor microenvironment and immune response analysis
The infiltration level of immune cell score and stromal cell score was obtained by using the "estimate" R package. The immune score and stromal score associated with the risk score or prognostic model were analyzed by Spearman correlation. Stem-cell-like features with risk score or prognostic model were measured by tumor stemness features extracted from the TCGA tumor samples. Six immune subtypes were defined to measure immune infiltrates. In tumor immune response, differences in prognostic model risk scores between the immune infiltrate subtypes obtained from the TCGA-LIHC were calculated by ANOVA analysis.

Chemotherapy sensitivity analysis
The NCI-60 database was accessed using the CellMiner interface (https://discover.nci.nih.gov/cellminer/) containing 60 different cancer cell lines from 9 different tumor types. Pearson correlation was used to analyze the association between the prognostic gene expression and drug sensitivity by using 263 drugs approved by the FDA or obtained from clinical trials (Supplementary Table 2).

Verification of mRNA expression of prognostic genes between HCC and adjacent non-tumorous tissues by qRT-PCR
The qRT-PCR experiments were performed to validate the mRNA expression levels of the five prognostic AGING genes in 20 paired HCC and adjacent non-tumorous tissues recruited from the First Affiliated Hospital of Wenzhou Medical University (Wenzhou, China). This study was approved by the Review of Ethics Committee in Clinical Research of the First Affiliated Hospital of Wenzhou Medical University. Total RNA from the HCC and adjacent normal liver tissue specimens was prepared with the Trizol reagent following the manufacturer's protocol (Servicebio). Using the RevertAid First Strand cDNA Synthesis Kit (Thermo), RNA was reverse-transcribed into cDNA. Gene expression was normalized to GAPDH. RT-PCR analysis was quantitated with FastStart Universal SYBR Green Master (Roche) by ABI StepOne (Applied Biosystems). The primer sequences were detailed in Supplementary Table 3. Each RNA sampling was performed in triplicate. To compare the expression levels between different samples, the relative expression of focal adhesion related genes was detected by the 2 −ΔΔCt method.

Verification of the protein expression of prognostic genes between HCC and adjacent non-tumorous tissues by IHC
IHC experiments were performed to validate the protein expression levels of the five prognostic genes in 10 paired HCC and adjacent non-tumorous tissues. All specimens were fixed with 10% formalin at room temperature, embedded in paraffin, and consecutively sliced into 4 µm sections. In brief, the tissue slices were firstly deparaffinized, and boiled in 10 mmol/L citrate buffer (pH = 6.4) for 10 min to retrieve the antigen. Then, the sections were treated in methanol containing 3% hydrogen peroxide to inactivate the endogenous peroxidase, and then with citrate buffer (pH = 6.0) to optimize antigen retrieval. The sections were incubated with 1% bovine serum albumin (BSA) in phosphatebuffered saline (PBS) for 30 min to block the unspecific binding. Besides, primary antibodies and HRPconjugated secondary antibodies were respectively applied to stain the slides. The detail information of antibodies is provided in Supplementary Table 5. After that, the slices were stained diaminobenzidine and counter-stained with hematoxylin. Finally, the sample was sealed, observed, and photographed.

Statistical analysis
Gene expressions between tumor and adjacent nontumorous tissues were compared by Wilcoxon test. Differences in proportions were compared by the Chisquared test. OS in different groups was detected by log-rank test and compared by Kaplan-Meier analysis. Independent predictors of OS were screened by both univariate and multivariate Cox regression analyses.
The ssGSEA scores of immune cells or immune pathway activities between high-and low-risk groups were compared by Mann-Whitney test. Spearman or Pearson correlation was used to explore the correlation between prognostic model risk scores or prognostic gene expression levels and stemness score, stromal score, immune score and drug sensitivity. Plots were created using R software (Version 3.6.3) with packages venn, igraph, ggplot2, pheatmap, ggpubr, corrplot, and survminer where appropriate. A twotailed p-value of <0.05 was considered to be statistically significant.

CONFLICTS OF INTEREST
The authors have declared that no conflicts of interest exists.