CELSR3 mRNA expression is increased in hepatocellular carcinoma and indicates poor prognosis

Data mining and collection

HCC patients in the TCGA and GEO cohorts meeting the following criteria were included in the study: (1) primary HCC samples (including hepatocellular carcinoma, hepatocholangiocarcinoma (mixed) and fibrolamellar carcinoma); (2) complete RNA-seq data. The exclusion criteria was samples that did not contain enough data for analysis.The gene expression data (423 cases, Workflow Type: HTSeq-Counts) and corresponding clinical information were obtained from The Cancer Genome Atlas (TCGA) Liver Hepatocellular Carcinoma (TCGA-LIHC) study of the official TCGA website (https://cancergenome.nih.gov/).The complete clinical data of the corresponding patients were obtained from cBioPortal (http://www.cbioportal.org/). A total of 373 HCC tissue specimens and 50 adjacent nontumor tissue specimens of HCC were included in the study. The Homo_sapiens.GRCh38.84.chr.gtf.gz file was downloaded from the Ensembl website (https://asia.ensembl.org/index.html), and ID conversion was performed using Perl (version 5.26.1). TCGA original HTSeq-Count data were processed using the Trimmed Mean of M-values (TMM) method for homogenization with edgeR. When an RNA had duplicate data, the average RNA expression was used. Microarray data were downloaded from the Gene Expression Omnibus (GEO) database. The GSE54236 original Series Matrix dataset with survival data was based on the GPL6480 (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F, Agilent Technologies, CA, Palo Alto) (submission date: Jan 21, 2014). The GSE60502 (submission date: Aug 18, 2014) datasets were based on GPL570 (Affymetrix Human Genome U133A Array Affymetrix, Inc, CA,Santa Clara). GSE41804 (submission date: Oct 24, 2012), GSE45436 (submission date: Mar 22, 2013), GSE62232 (submission date: Oct 09, 2014), and GSE6764 (submission date: Jan 17, 2007) datasets, which were based on GPL570 data (Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix, Inc, CA, Santa Clara). For the CEL expression profiles, Robust Multiarray Average (RMA) normalization was performed using the affy package. The expression level of CELSR3 mRNA was converted using the log2 value for further analysis. The raw data in the dataset were annotated to obtain the gene expression levels and the average expression values of probes were considered as the expression values of the corresponding genes. We combined the HCC microarray GSE60502 and GSE62232 in the GEO database using Perl (version 5.26.1) and performed batch normalization using the sva package of R. After combination and normalization, 99 HCC patients were divided into high and low expression groups according to the median value of CELSR3 expression for GSEA detection and verification. RNA expression (RNA-Seq) data for CELSR3 in established HCC cell lines were accessed on 01∕02∕19 (n = 25) from the Cancer Cell Line Encyclopedia (CCLE) database (https://portals.broadinstitute.org/ccle/about); RNA expression values were reported in reads per kilobase of transcript per million mapped reads (RPKM). We divided the expression of CELSR3 in the HCC cell lines into high and low expression groups using the median value, and then performed GSEA verification to observe the enrichment of CELSR3 differential expression in the HCC cell lines.The GSEA results were presented in the form of multipleGSEA using the Plyr, grid, gridExtra and ggplot2 packages of R.

Gene set enrichment analysis

GSEA is a computational method that determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states (Subramanian et al., 2005). In this study, HCC samples from TCGA data were divided into the CELSR3 high mRNA expression group and the low CELSR3 mRNA expression group based on the median value of CELSR3 mRNA expression. GSEA3.0 was adopted for GSEA. The functional gene set files “c2.cp.biocarta.v6.2.symbols.gmt”, “h.all.v6.2.symbols.gmt” and “c5.all.v6.2.symbols.gmt” were used to summarize and elucidate specific and well-defined biological states or processes. The number of substitutions per analysis was set at 1,000, and gene sets with P < 0.05 and a false discovery rate (FDR) < 0.25 were recognized as a significantly enriched.

Statistical analysis

SPSS statistical software, version 19.0 (SPSS Inc., Chicago, IL, USA) and STATA 12.0 (Stata Corporation, College Station, TX, USA) were used for statistical analysis. The ggplot2 and pROC packages in the statistical software R (R Core Team, 2018) were used for graph generation. Discrete variables are expressed using a box plot to measure expression differences. The chi-square test was used to analyze the relationship between CELSR3 mRNA expression and clinical data. The Kaplan–Meier curve showed that clinicopathologic characteristics were associated with overall and relapse-free survival. Univariate Cox analysis was used to select relevant variables, and subsequently, multivariate Cox analysis was used for prognostic analysis of CELSR3 mRNA expression with regard to the overall and relapse-free survival (RFS) rate of HCC patients. The cutoff value was determined by the median value of CELSR3 mRNA expression. P < 0.05 was considered statistically significant.

Patient characteristics

The clinical and gene expression data for the 373 cases of primary HCC were downloaded from TCGA. As shown in Table 1, the median age at diagnosis was 61 years. In our study cohort, after excluding data that is not available, 14.9% of the tumors were well differentiated, 48.4% were moderately differentiated, and 36.7% were poorly or not differentiated. The tumor status included 234 cases without tumors (67.6%) and 112 cases with tumors (32.4%). There were 172 cases (49.3%) of stage I primary HCC, 87 cases (24.9%) of stage II HCC, 85 cases (24.4%) of stage III HCC, and 5 cases (1.4%) of stage IV HCC. There were 326 cases without and 18 with residual tumors. Of the subjects who were alive at the last follow-up, the median survival time was 19.42 months (range: 0-120.73 months).

High CELSR3 mRNA expression in HCC

HCC tissues and normal control tissues in TCGA data were compared using a box plot. As shown in Fig. 1A, the expression of CELSR3 mRNA was higher in HCC tissues (P = 0.000). Moreover, there were also different CELSR3 expression levels in the groups classified according to age (Fig. 1B, P = 0.023), tumor status (Fig. 1C, P = 0.030), stage (Fig. 1D, P = 0.023), and vital status (Fig. 1E, P = 0.027). The results indicated that CELSR3 is highly expressed in HCC tissues, and the expression of CELSR3 was high in older patients, in those with advanced tumor status or TNM stage III, and in those who had died. It is interesting to note that the expression of CELSR3 decreases in stage IV. When the GEO data were validated, CELSR3 mRNA expression differed for different types of liver tissues in the GSE41804 (Fig. 1F, P = 0.000), GSE45436 (Fig. 1G, P = 0.000), and GSE6764 (Fig. 1H, P = 0.000) datasets. It is worth noting that the expression of CELSR3 mRNA in liver cirrhosis tissues was lower than in normal and HCC tissues.

Diagnostic function of CELSR3 mRNA

To verify the diagnostic value of CELSR3 in HCC, the ROC curve was used to analyze the AUC of CELSR3 expression associated with clinical pathological parameters of different HCC patients. As shown in Fig. 2A, CELSR3 could effectively distinguish normal liver tissue from HCC tissue. In the receiver operating characteristic (ROC) curve analysis of CELSR3, the area under the curve (AUC) was 0.954, indicating high diagnostic and predictive values of CELSR3. Furthermore, we performed ROC curve analysis of CELSR3 expression in different TNM stage subgroups of HCC patients; our results indicate that high expression of CELSR3 may be a promising diagnostic marker. In the subgroup analysis of stage, the diagnostic values at different stages were also high (AUC: stage I = 0.946; stage II = 0.964; stage III = 0.967; stage IV = 0.936; Figs. 2B–2E).

Relationship between CELSR3 mRNA expression and clinical features of HCC

To validate the relationship between CELSR3 expression and clinical pathological features in HCC patients, we further analyzed the expression level of CELSR3 in HCC patients at different clinical stages. As shown in Table 2, the relationship between CELSR3 mRNA expression and clinical features indicated that CELSR3 mRNA expression was significantly correlated with age (P = 0.025), tumor status (P = 0.022), clinical stage (P = 0.003), T classification (P = 0.010), vital status (P = 0.001), and relapse (P = 0.005). Collectively, these data indicate that the expression of CELSR3 is associated with various important clinical pathological features of HCC.

Survival results and multivariate analysis

The Kaplan–Meier survival curve and log-rank test were used to evaluate the relationship between CELSR3 mRNA expression and overall (Fig. 3) and relapse-free survival (Fig. 4). All HCC patients were divided into either the CELSR3 high expression group or CELSR3 low expression group based on the median value of CELSR3 mRNA. We analyzed the OS time of 372 patients. The result showed that the overall survival was poor in patients with high CELSR3 mRNA expression (Fig. 3A; P = 0.000). Unexpectedly, RFS analysis produced similar results (Fig. 4A; P = 0.000). These results indicated that high CELSR3 expression in patients with HCC can predict a poor prognosis. To further confirm the prognostic value of CELSR3 in patients with HCC, we performed OS and RFS analysis in subgroups of HCC patients. The subgroup analysis indicated that the overall survival was poor in patients with high CELSR3 mRNA expression and an age <60 years or ≥ 60 years, male sex, histological grade G1/G2, stage I/II, N0 stage, M0 stage, and R0 stage (Figs. 3B–3I), and relapse-free survival was poor in the group with high CELSR3 mRNA expression and an age <60 years, female sex or male sex, histological grades G1/G2, stage I/II, T1 stage, N1Nx stage, M0 stage, and R0 stage (Figs. 4B–4J). The survival analysis validated by the GSE54236 data is shown in Fig. 3J (P = 0.007). These results demonstrated that high expression of CELSR3 can function as a prognostic biomarker of OS and RFS in subgroups of HCC patients. In different clinical subgroups of HCC patients, the prognostic value of CELSR3 varies, which can guide our clinical practice and deserves further study.

Univariate analysis showed that high expression of CELSR3 mRNA was significantly correlated with poor overall survival and other variables associated with a reduced overall survival rate, including tumor status, stage, T stage, N stage, M stage, and residual tumor. Multivariate analysis using the Cox proportional hazards model indicated that high expression of CELSR3 mRNA (HR = 1.88, P = 0.004) and residual tumor (HR = 1.36, P = 0.033) were independent prognostic factors for the overall survival of HCC patients (Table 3).

Identification of CELSR3-related signal transduction pathways by GSEA using TCGA cohort

To identify the differentially activated signaling pathways in HCC, gene expression enrichment analysis was performed between datasets with low or high CELSR3 mRNA expression. GSEA revealed significant differences (FDR < 0.25, NOM P-value <0.05) in the enrichment of the MSigDB Collection (c2.cp.biocarta.v6.2.symbols.gmt, h.all.v6.2.symbols.gmt, and c5.all.v6.2.symbols.gmt). Based on the normalized enrichment scores (NESs), the most significantly enriched signal transduction pathways were selected (Figs. 5A–5C; Tables 4 and 5). Figures 5A–5C shows that the cell cycle pathway, MCM pathway, ATR and BRCA pathway, ATM pathway, biocarta G1 pathway, E2F targets, hallmark G2/M checkpoint, and spermatogenesis, etc., were differentially enriched in phenotypes with high CELSR3 expression.

Validation using meta-GEO HCC cohort from GEO database and cell line data from CCLE

A total of 99 HCC specimens were obtained after combining the GSE60502 and GSE62232 microarray and batch normalization treatment. We performed enrichment analysis of meta-GEO HCC cohort data and cell line data. The results showed that after choosing the c2.cp.biocarta.v6.2.symbols.gmt gene set, two common enriched pathways,

biocarta G1 pathway, ATR and BRCA pathway, were present in both the meta-GEO cohort and the TCGA cohort of the high CELSR3 mRNA expression phenotype (Fig. 5D; Table 4). After choosing the h.all.v6.2.symbols.gmt gene set, there were two common enriched pathways, hallmark G2/M checkpoint and hallmark spermatogenesis, in the meta-GEO cohort and the TCGA cohort of the high CELSR3 mRNA expression phenotype (Fig. 5E; Table 4), and one common enriched pathway, hallmark DNA repair, in the meta-GEO cohort and the HCC cell line cohort of the high CELSR3 mRNA expression phenotype (Fig. 5G; Table 4). After choosing the c5.all.v6.2.symbols.gmt gene set, four pathways, GO somatic diversification of immune receptors, GO endonuclease activity, GO DNA repair complex and GO somatic cell DNA recombination, were differentially enriched in the meta-GEO cohort, the HCC cell line cohort and the TCGA cohort of the high CELSR3 mRNA expression phenotype (Figs. 5C, 5F, 5H; Table 5).