Screening and Identification of Potential Biomarkers for Hepatocellular Carcinoma: An Analysis of TCGA Database and Clinical Validation

Introduction

HCC is one of the most common malignancies worldwide with an increase in incidence.^1,2 Currently, the most widely used diagnostic methods for HCC screening include biannual abdominal ultrasound examination and serum alpha-fetoprotein measurement.³ In the United States, all focal liver lesions suspected during ultrasound examination are verified using computed tomography and/or magnetic resonance imaging.⁴ The application of these methods makes the diagnosis of HCC a lot more accurate, in which the sensitivity of diagnosis is up to 89% and the specificity reaches as high as 99%.⁵ Unfortunately, if the lesions are <1 cm in diameter, the outcome of diagnosis may not be as precise (could be as low as 34%).⁶ Thus, it is of great significance to identify novel biomarkers of HCC that provides more opportunities for diagnostic indicators or gene-targeted therapy of HCC.

With the application of high-throughput sequencing technology, a large amount of cancer data has been collected and is now available to the medical research community.^7,8 By extracting and analyzing these valuable data, it is possible to identify potential biomarkers that can improve the detection of cancer.⁹ RF is an ensemble algorithm that combines many weak predictors with good local performance to yield an ensemble predictor with strong global performance.^10,11 In cancer research, RF could be used to predict the risks of cancer and the susceptibility of patients. And the feature importance of RF can be used to identify the significance of HCC tumor-related genes.

In this study, we established an effective classifier using the RF algorithm and thereby analyzed HCC transcriptome data plus corresponding control data obtained from TCGA database. Consequently, we could identify several potential biomarkers for HCC, which were verified in cancer cell lines and situ tumor samples. The survival curve revealed a negative correlation between the potential biomarkers and survival rates of patients.

Materials and Methods

Data Collection

HCC RNA-Seq data were downloaded from the TCGA database using The GDC Data Portal (https://portal.gdc.cancer.gov/). A total of 23,216 feature genes were revealed. Overall, 373 HCC tissue samples and 50 adjacent non-tumor tissue samples were included in the study.

Identification of Potential Biomarkers Using the RF Algorithm

RF is an ensemble algorithm that is composed of many decision trees.¹² To classify a new instance, each decision tree provides a classification for input data and RF collects the classifications and chooses the most voted prediction as to the result, which has high accuracy values and high prediction power.¹³ An RF classifier was constructed based on all genes from the TCGA HCC database to identify potential biomarkers for HCC. All samples were stochastically rearranged and divided into 10 parts: seven of the 10 parts were used as the training set to acquire the threshold parameter for a training classifier, while the other parts were used as the test set. Subsequently, the training classifier was used to predict the accuracy within the test set by calculating the false-positive rate and false-negative rate. Receiver operating characteristic curve analysis was performed to evaluate the classification performance and robustness of the prediction classifier, which represented the ability of a feature to contribute to the prediction of target variables in the model, was calculated by a total decrease in tree-node impurities from splitting on the predictor feature variable and averaged over all sub-trees in the RF classifier. The importance measures of RF were often used for feature selection for high-dimensional data settings, which, in particular, makes RF attractive for bioinformatics and related fields wherein identifying a subset of relevant predictors from a large set of candidate predictors is a major challenge (ie, the “small n large p” problem). Features with FI of >5% were regarded as candidate genes that could be potential biomarkers. (The specific code was shown in Code S1, S2, and S3.)

Three-Level PPI Network Construction

The candidate genes selected from the RF classifier were mapped to the STRING database (http://www.string-db.org/) to find their first- and second-order outcome-associated proteins.¹⁴ PPI data were downloaded from the STRING database and visualized by using Cytoscape 3.6.1. Network topology characteristics were calculated using the Network Analyzer plugin.¹⁵

GO Functional and KEGG Pathway Enrichment Analysis

We performed GO and KEGG pathway enrichment analyses to determine potential biological functions and pathways of the protein mentioned in a three-level PPI network using the Database for Annotation Visualization and Integrated Discovery online tool (https://david.ncifcrf.gov/). p<0.05 was considered statistically significant.

Survival Analysis

Survival rates were assessed using the Kaplan–Meier method and survival curves were plotted using an online tool (http://gepia.cancer-pku.cn/detail.php). Comparisons between study groups were made using a Log-rank test; log-rank p<0.05 was considered significant.

Cells

The hepatocellular carcinoma HepG2 cells and normal liver LO2 cells were purchased from the American Type Culture Collection (http://www.atcc.org). The human HCC cell lines (BEL-7402 and SMMC-7721) were bought from the Shanghai Institute of Biotechnology and Cell Biology (Shanghai, China). All cell lines were maintained in a humidified 5% CO₂ incubator in low-glucose Dulbecco’s Modified Eagle’s Medium (DMEM) culture medium containing 10% fetal bovine serum (FBS), 100 U/mL penicillin, and 100 μg/mL streptomycin. For sub-culturing, cells were detached using 0.05% trypsin–EDTA treatment at 37°C. And all cell lines were regularly checked by Quick Cell Mycoplasma Rapid Detection Kit (Shanghai Life iLab Biotech Co., Ltd, AC16L061) to ensure they are authentic and keep in good shape.

Tissue Specimens

Five paired fresh HCC specimens and adjacent non-tumor liver tissues were obtained from patients who underwent HCC resection at the Department of Abdominal Surgery, Cancer Center of Guangzhou Medical University (Guangzhou, China) in 2019. All tissues were used for Western blotting analysis. This study was approved by the Ethics Committee of the Cancer Center of Guangzhou Medical University. The study was conducted in accordance with Declaration of Helsinki. Written informed consent was obtained from every patient.

Western Blotting

Cells (LO2, BEL-7402, SMMC-7721, and HepG2) and tissues (HCC and adjacent non-tumor tissues) were lysed in a loading buffer (1×). After lysis, protein extracts (10 μL/lane) were subjected to 12% SDS-PAGE and transferred to polyvinylidene fluoride membranes. The membrane was blocked in TBST (0.05% Tween 20 in Tris-buffered saline) containing 5% milk for 2 h and then probed with antibodies for NOX4 (Absin Bioscience, Shanghai, China) and FLVCR1 (ABclonal Technology, Wuhan, China), followed by incubation at 4°C for overnight. The membrane was then washed three times with TBST, followed by incubation with secondary antibodies for 1 h at room temperature. Immunoreactive bands thus obtained were visualized by enhanced chemiluminescence.

Immunofluorescence Assay

Cells (1 × 10⁵ cells/well) were seeded onto a copolymer culture dish. The cells were fixed with 4% paraformaldehyde for 20 min at room temperature. For permeabilization, the cells were incubated with a blocking solution (9 mL PBS, 1 mL FBS, and 30 μL Triton X-100) for 3 h at room temperature. Subsequently, they were incubated with primary antibodies in a blocking solution at 4°C overnight, and then incubated with fluorescent secondary antibodies (1/400) at room temperature for 2 h; finally, they were stained with Hoechst 33342 for 10 min. Between each step, the cells were washed five times with PBS. Images were captured using the LSM 800 confocal microscope platform. Images were obtained with a laser scanning confocal microscope (Karl Zeiss, Jena, Germany).

IHC

Twenty-four cases of HCC and adjacent tissue chips were purchased from Wuhan Iwill Biological Technology Co., Ltd (Wuhan, China). Immunohistochemical determination for NOX4 (dilution 1:200) and FLVCR1 (dilution 1:400) was performed according to standard procedures using an antigen retrieval technique in citrate buffer. Follow-up IHC was performed with a ready-to-use high-efficiency immunohistochemical secondary antibody kit (Absin Bioscience, Shanghai, China), according to manufacturer instructions. Finally, we determined protein expression by optical microscopy. Two observers performed independent examinations to record histopathological features and specific immunoreactivity (IR). Staining was semi-quantitatively assessed by assigning scores to the intensity of IR and the proportion of positively stained cells. IR intensity was divided into four categories: 0 (no IR); 1 (weak IR); 2 (moderate IR); and 3 (strong IR). The proportion of positively stained cells was divided into four groups: 0 (no tumor cells showing IR); 0.33 (<1/3 of the tumor cells showing IR); 0.67 (1/3–2/3 of the tumor cells showing IR); and 1 (>2/3 of the tumor cells showing IR). The product of the two values represented the overall IR score.¹⁶

Results

Top 10 Genes with the Best Feature Importance Were Screened as Candidate Genes

To identify potential cancer-related genes for HCC, an RF classifier was established (Figure 1A). The performance of the RF algorithm was compared against that of the decision tree and gradient boosting decision tree algorithms. Our experimental results demonstrated that the RF classifier relatively exhibited the best performance, with an area under the curve (AUC) of 0.9974 and a comprehensive evaluation index (F1) score of 99.12% (Figure 1B, Table 1). Then, we identified 10 genes (NOX4, PAMR1, NAT2, STIL, NUDT10, C16ORF59, SFTA1P, FLVCR1, ZADH2, and VMO1) with high feature importance (FI > 5%) as candidate genes. The importance values of each candidate gene are shown in Figure 1D. Subsequently, we confirmed the importance of the 10 candidate genes by using them alone as the basis of classification (features) in order to classify all samples and to make a prediction. As a result, the 10 genes showed high classification performance with AUC of 0.9981, which was even higher than the AUC value of the original RF classifier (Figure 1C).

Table 1 Classification Results of Three Algorithms

Figure 1 Ten candidate genes of HCC identified using RF algorithm. (A) Study workflow. (B) The receiver operating characteristic curve threshold was used for measuring algorithm quality. (C) The receiver operating characteristic curve of 10 important candidate genes under the RF classifier. (D) Importance value of top 10 genes. Abbreviations: FI, feature importance; AUC, area under the curve; PPI, protein–protein interaction; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; IHC, Immunohistochemistry.

Candidate Genes are Involved in Cancer-Related Pathways

Since the expression levels of the 10 candidate genes significantly different between HCC tissue samples and adjacent normal tissue samples, the biological roles of these genes warranted further investigations (Figure 2B). We considered the proteins encoded by all candidate genes as driver proteins and constructed a three-level PPI network. As one of the candidate genes, named SFTA1P, is a non-coding RNA, the first PPI network level only contained nine candidate proteins. Besides, the three-level PPI network also had 1280 outcome-associated proteins that linked with the nine candidate proteins, included 161 first-order outcome-associated proteins in the second PPI network level and 1119 second-order outcome-associated proteins in the third PPI network level (Figure 2A). So there were 1289 outcome-associated proteins in the three-level PPI network. In GO analysis, these 1289 proteins were mainly involved in mitotic nuclear division, nuclear division, and organelle fission at the biological process level. Besides, at the cellular component level, chromosomal region, centromeric region, and condensed chromosome were the three most significant categories. And at the molecular function level, peptide receptor activity, tubulin binding, and G protein-coupled peptide receptor activity were the most significantly enriched (Figure 2D). Furthermore, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using the three-level PPI network proteins. As a result, the proteins were significantly enriched in cancer-associated pathways, which meant that the nine candidate proteins participated in cancer progress (Figure 2C).

Figure 2 Association of the 10 candidate genes with tumor development. (A) Three-level protein–protein interaction network results. Red boxes represent the candidate genes, and yellow and blue boxes represent first- and second-order outcome-associated genes, respectively. The width of edges varies with their confidence. (B) A heat map depicting different expression levels of the 10 candidate genes in hepatocellular carcinoma tissues versus adjacent non-tumor tissues. Green represents low expression and red represents high expression. (C) KEGG pathway analysis. (D) GO enrichment analysis results. p<0.05 indicated statistical significance.

Upregulated Expression Levels of NOX4 and FLVCR1 in HCC

As reported above, we predicted the possible upregulation of 10 genes in HCC tissue samples compared with adjacent normal tissue samples. To validate our prediction, we selected two representative genes, ie, NOX4 and FLVCR1, the expression of which has rarely been reported in liver (LO2) and HCC (BEL-7402, SMMC-7721, and HepG2) cells. In comparison with LO2 cells, HCC cells showed a marked increase in the protein expression levels of NOX4 and FLVCR1 (Figure 3A). Immunofluorescence assay was conducted to confirm the Western blotting results (Figure 3B). Immunohistochemistry (IHC) also revealed high expression levels of NOX4 and FLVCR1 in HCC tissues but weak expression levels in adjacent non-tumor tissues (Figure 4A). The expression levels of NOX4 and FLVCR1 in fresh tissue samples were also examined by Western blotting, which further confirmed that NOX4 and FLVCR1 expression levels in HCC tissues were higher (Figure 4B).

Figure 3 Upregulation of NOX4 and FLVCR1 in HCC cells. (A) Expression levels of NOX4 and FLVCR1 in normal liver cells (LO2) and HCC cells (BEL-7402, SMMC-7721, and HepG2), as determined by Western blotting. ***p<0.001 compared with LO2 cells. (B) Immunofluorescence assay of NOX4 and FLVCR1. NOX4 and FLVCR1 are labeled in red; the nucleus was stained with Hoechst (blue).

Figure 4 Upregulation of NOX4 and FLVCR1 in HCC tissues. (A) Expression levels of NOX4 and FLVCR1 in HCC tissues and adjacent non-tumor tissues, as evaluated using immunohistochemical staining (×100). Expression levels of both NOX4 and FLVCR1 were higher in HCC tissues than in adjacent non-tumor tissues. (B) Western blotting showing high expression levels of NOX4 and FLVCR1 in HCC tissues. *p<0.05, **p<0.01, ***p<0.001 (data compared to normal tissues). Number # indicated liver samples from different HCC patients. Abbreviations: T, tumor tissues; N, non-tumor tissues.

High Expression Levels of NOX4 and FLVCR1 Associated with Lower Survival Rates

Next, we analyzed the clinical significance of NOX4 and FLVCR1 in patients with HCC. Kaplan–Meier analysis indicated patients with HCC that highly expressed NOX4 or FLVCR1 had relatively lower survival rates (a hazard ratio >1.5 for these patients; Figure 5), which implied these two genes, in addition to serving as biomarkers for HCC cells, may also be involved in tumor development and progression.

Figure 5 Survival analysis of NOX4 and FLVCR1. Kaplan–Meier analysis of the relationship between overall survival and expression levels of NOX4 (A) and FLVCR1 (B); Log-rank p< 0.05 indicated statistical significance.

Discussion

Considering the increasing speed and reduced cost of the sequencing technology, a huge amount of transcriptomic data has been generated, facilitating the comprehensive elucidation of gene expression profiles at the transcriptional level in cancer cells. Researchers working on cancer-driven topics can now conveniently acquire transcriptomic data from TCGA to perform secondary analysis of gene expression via bioinformatic tools. However, the traditional approach, eg, analyzing differentially expressed genes, generally highlights non-specific biomarkers considering the involvement of complex interactions (Figure S1, Table S1). Herein we initially constructed a classifier to analyze the HCC RNA-Seq data obtained from the TCGA database and identified 10 genes with high significance in distinguishing HCC tissues from adjacent non-tumor tissues. The upregulation of NOX4 and FLVCR1 in HCC cell lines or clinical tumor tissues confirmed our prediction, illustrating that the RF classifier had an excellent performance. PPI network analyses indicated that the predicted genes were involved in complex tumor growth-related pathways, suggesting that they had a strong correlation with HCC development. Consistently, Kaplan–Meier analysis illustrated that patients with HCC with high expression levels of NOX4 or FLVCR1 had a lower survival rate, suggesting that the upregulation of these genes represents a poor prognosis for these patients. NOX4 and FLVCR1 may thus be involved in tumor development and progression in addition to serving as biomarkers for HCC cells.

We identified 10 candidate genes using the RF classifier through the integrated analysis of the TCGA database. Next, we determined the clinical significance of the 10 candidate genes in patients with HCC. As evident from the survival analysis, NOX4 and FLVCR1 showed better performance as biomarkers than the other eight candidate genes (Figure S2). Thus, we selected NOX4 and FLVCR1 to further examine their expression levels in HCC cells. Both HCC cell lines and tumor tissues showed the intracellular overproduction of NOX4 and FLVCR1, presenting them as potential biomarkers. To determine whether the predicted genes correlated with tumor development and progression, we constructed a three-level PPI network using the candidate genes as driver proteins that included first- and second-order outcome-associated proteins (Figure 2A). The vector illustration was shown in Figure S3. KEGG analysis showed that the proteins were involved in complex cell signal pathways associated with tumor growth and metastasis.^19–21 This result suggested that the candidate genes can not only be regarded as biomarkers for HCC but also act as tumor-promoting factors, which may be potential therapeutic targets for HCC.

Consistent with our hypothesis, NOX4 as an oxidoreductase has been reported to catalyze ROS production by NADPH and promote prostate cancer and glioblastoma.^22,23 Moreover, Roy et al elucidated the potential mechanisms of NOX4-related damage that contribute to the progression of solid malignancies; they reported that the over-expression of functional NOX4 in numerous tissues helps explain tissue injury and DNA damage due to ROS that accompanies pre-malignant conditions.²⁴ Although FLVCR1 encoding a ubiquitous heme exporter is reportedly associated with leukemia,^25,26 the relationship between FLVCR1 and HCC has not been reported as yet. In synovial sarcoma, FLVCR1 was demonstrated to promote the proliferation and tumorigenicity through inhibiting apoptosis and autophagy.²⁷ Similarly, as apparent from the results of our study, FLVCR1 may be an HCC tumor-promoting factor, which reveals a potential diagnostic and therapeutic target for HCC.

Conclusions

Herein we identify two potential biomarkers for HCC, namely NOX4 and FLVCR1, which may unravel novel diagnostic and therapeutic opportunities.