This study was approved by the Institutional Review Board of The First Affiliated Hospital. All tissue samples used in the current study were obtained from patients with HCC who underwent a partial hepatectomy at the First Affiliated Hospital, Zhejiang University School of Medicine. Each specimen was reviewed by a board-certified pathologist to confirm that the frozen section was histologically consistent with tumor or non-tumor tissues. Written informed consent was obtained from each patient.

Paired tumor and adjacent non-tumor tissue samples from 33 HCC patients were subjected to WGBS on the Illumina X Ten platform with the procedures described in our previous WGBS paper (29). Briefly, a 200-ng genomic DNA sample was sheared to about 300-bp fragments by sonication. Then, DNA fragments were subjected to end-repair, addition of adenosine to the 3′ end, and TruSeq adaptor ligation (Illumina, San Diego, CA, USA). Bisulfite conversion was implemented via the EZ DNA Methylation kit (Zymo Research, Irvine, CA, USA) according to the manufacturer’s protocol. After that, bisulfite-converted DNA was enriched through several cycles of PCR amplification using the KAPA HiFi HotStart uracil DNA polymerases (Kapa Biosystems, Boston, MA, USA). The PCR conditions were set as 45 s at 98°C followed by 10 cycles at 98°C for 15 s, 65°C for 30 s, 72°C for 30 s, ending with 72°C for 1 min. The quality of each WGBS library was assessed by Qubit 2.0 (Life Tech, Carlsbad, CA, USA) and an Agilent 2100 Bioanalyzer. Finally, 150-bp pair-end sequencing was conducted on the Illumina X Ten sequencing platform. High-throughput mRNA-seq was performed for each WGBS sample. Similarly, all 66 RNA samples with high quality (RIN ≥ 7) were applied to the Illumina X Ten platform for sequencing. Specifically, total RNA was extracted and purified using the RNeasy Micro Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. The quality of RNA was assessed via an Agilent 2100 Bioanalyzer. Libraries for poly(A)+ RNA were prepared according to the Illumina standard protocol. Constructed libraries were sequenced on HiSeq X Ten platform by WuXi AppTec (Wuxi, Jiangsu China).

For raw reads from RNA-seq, Cutadapter (30) (v.1.12) and Trimmomatic (31) (v. 0.33) were applied for adapter removal and trimming of low-quality sequences, followed by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) for a quality check. After quality control, clean reads were submitted as the input of Kallisto (32) (v.0.44) for abundance quantification of transcripts based on a gene model download from the GENCODE (v. 29) (33). Normalized expression (transcripts per million [TPM] reads) of each gene was summarized from the transcript level via the R package tximport (v. 1.6.0) (34). Non-expressed and low-expressed genes (defined as those with TPM < 0.01 among more than half of the total 66 samples) were excluded from downstream analysis. Clean reads were also aligned against the reference genome by STAR (v. 2.5.2a) (35), and the resulting BAM files were utilized as input for enhancer RNA expression quantification via bedtools (v.2.27.1) (36). Besides, the gene count outputs of STAR were used as inputs of DESeq2 (v.1.18.1) (37) for identification of differential expressed genes (FDR < 5% and |Log₂FoldChange| > 0.5). The clean WGBS reads that passed preprocessing were aligned with the hg38 reference genome using Bismark (v. 0.16.1) (38) with default parameters. The harvested count information from each strand was then combined. As recommended by the R package DSS (v.2.26.0) developer (39), the smoothing approach was adopted for estimation of smoothed methylation level for all 28.9 million CpGs with default parameters.

In order to identify overall significant differential methylation between all tumor and non-tumor samples, a combined Baumgartner–Weiβ–Schindler (BWS) test (40) was applied to carry out age-adjusted differentially methylated loci (DML) detection via the R package BWStest (v.0.2.2). We divided all the 33 HCC patients into three age groups: “young” (age < 55 years; n = 10), “medium” (55 < age ≤ 65; n = 13), and “old” (age > 65; n = 10). A single BWS test was performed for each age group on every CpG to obtain two individual BWS p-values (p _left and p _right). Afterward, three one-sided p-values were combined as statistic T_left (or T_right) = -2* Σlog₁₀(p _left [or p _right]), and a new statistic T was defined as max(T_left, T_right) (40). The empirical distribution of the T statistics of combined BWS test was determined by 2.0 × 10⁸ time permutations. At last, an overall empirical p-value was estimated as the combined BWS p-values for each CpG. CpG with a combined BWS p-value < 1.0 × 10^-5 was identified as DML for subsequent differentially methylated loci (DMR) calling.

Tumor-associated DMRs were determined by R script with the following two steps: 1) DML were combined into pre-DMRs if the distance between neighbor CpGs was < 200 bp; and 2) all CpGs located between the start and the end of each pre-DMR were included as a final DMR. The arithmetic mean of T statistics for all CpGs in each DMR was calculated for estimating the empirical combined BWS p-value for each DMR. Group-level methylation was estimated as the arithmetic average of DNA methylation of all CpGs in each corresponding DMR.

Annotation of the genomic location of each identified DML and DMR was realized by using the function “findOverlaps” in the R package GenomicRanges (v.1.30.3) (41). Specifically, location annotations were determined according to overlaps (>1 bp) between the range of each DML/DMR and all known genomic regions. These genomic region annotations were defined on the basis of the gene model downloaded from the GENCODE, which included promoter (upstream 1,500 bp and downstream 500 bp from the TSS of each transcript), exon, intron, 5′-UTR, 3′-UTR, and intergenic region.

Functional annotation was performed for all identified DMRs to search for potentially active promoter- and enhancer-like DMRs. Specifically, liver-active promoter/enhancers were obtained according to candidate regulatory element annotation of eight liver-relevant histone ChIP-seq samples, which included five tissue samples (i.e., two tumor and adjacent samples from two HCC patients and one normal liver sample) from a recent integrative epigenomic HCC study (42) and three liver-relevant samples (i.e., one adult liver tissue, one hepatocyte, and one HepG2 cell sample) from the ENCODE database (43). The regulatory element annotation of these five liver tissue samples was based on 10-state chromHMM annotations, whereas those three ENCODE samples were based on five-state candidate regulatory element annotations. Hence, active liver promoters were composed of regions annotated as “activeTSS” or “activePromoter” (refer to regions with both H3K4me3 and H3K27ac peaks) in any of those five tissue samples and “promoter-like cRE” (refer to regions with both H3K4me3 and DNase peaks) in any of those three ENCODE samples. Similarly, active liver enhancers included regions annotated as “activeEnhancer” (regions with both H3K4me1 and H3K27ac peaks) in any of those five tissue samples and “enhancer-like cRE” (regions with both H3K27ac and DNase peaks) in any of those three ENCODE samples. Finally, DMRs that overlapped with at least one active promoter/enhancer in any sample were identified as promoter-like or enhancer-like DMRs. Besides, promoter-like DMRs also included DMRs that were annotated as promoters in the genomic location annotation. For identified promoter/enhancer-like DMRs, their promoter/enhancer-like activity scores were calculated as the number of ChIP-seq liver samples in which they were annotated as promoter/enhancer-like regulatory elements. Moreover, as a supplement of annotated enhancers in the public domain, we estimated enhancer RNA (eRNA) expression, which was reported to be a reliable indicator of enhancer activity (44, 45), for all intergenic DMRs via bedtools (36) for identification of potential novel enhancers with active enhancer expression (count of reads ≥ 3 in at least one third of tumor or adjacent HCC samples) among our HCC samples.

To investigate the effect of the identified tumor-associated DMRs in transcriptional regulation, we performed integrative genomic analysis by integrating the paired DNA methylomic and transcriptomic data from our original HCC cohort and histone ChIP-seq data from the public domain. To ameliorate potential false-positive findings, we focused only on DMRs (|Δ_methylation| ≥ 0.15) in candidate active regulatory elements, i.e., active promoters and enhancers in the liver. As mentioned earlier, DMRs that overlapped with active promoters or enhancers were identified as promoter/enhancer-like DMRs. For promoter-like DMRs, nearby genes (TSS ≤ 2 kb away from the DMR start or end site) were tested for correlation between mRNA expression and methylation of DMR. Only differentially expressed genes (FDR < 5% and |Log₂ FoldChange| > 0.5) with a BH-corrected Spearman correlation p-value < 0.05 were designated promoter-like DMR-associated genes.

As for genic enhancer-like DMRs, nearby genes whose distances from DMR were < 100 kb were examined for Spearman correlation, and differentially expressed genes with a BH-corrected correlation p-value < 0.05 were identified as enhancer-like DMR-associated genes. Regarding those identified intergenic enhancer-like DMRs, neighbor genes within ± 0.5 Mb of those intergenic enhancer-like DMRs were screened for differentially expressed genes with significant DMR-eRNA–gene triple correlation (simultaneous significant Spearman correlation between DMR methylation and eRNA expression, between eRNA and gene expression, and between DMR methylation and gene expression). Additionally, enhancer-like DMR-associated target genes were ruled out of genes that were also identified as promoter-like DMR-associated genes. Those target genes associated with promoter-like DMRs, genic enhancer-like DMRs, or intergenic enhancer-like DMRs were defined as DMR-DEGs.

Genic- and intergenic enhancer-like DMR-associated genes were first combined as enhancer-like DMR-DEGs, along with promoter-like DMR-DEGs which were then used as inputs for enrichment analysis via Metascape (v.3.5) (46).

Strict screening procedures were imposed to identify DMR-DEGs with high confidence and a low possibility of false-positive results, which would be more applicable to downstream validation. For promoter-like DMRs and genic enhancer-like DMRs, we conducted another genomic location and functional annotation with a much stricter criterion. Specifically, only DMRs overlapping with at least 80% of a promoter region, an active promoter, or an active enhancer were defined as high-confidence promoter/enhancer-like DMRs. Afterward, all DEGs associated with those high-confidence promoter/enhancer-like DMRs were identified as high-confidence DMR-DEGs. As for intergenic enhancer-like DMRs, only target DEGs that were highly positively co-expressed (Spearman correlation coefficient ρ ≥ 0.7) with corresponding eRNAs passed the screening and were included as high-confidence DMR-DEGs.

In silico replication was conducted for each high-confidence DMR–DEG pair using methylomic and transcriptomic data of the TCGA-LIHC cohort. The clinical phenotypes, DNA methylation, and gene expression dataset of the TCGA-LIHC cohort were downloaded via the RTCGA (v.1.22.0) R package (47). The 450 k array methylation data were updated to hg38 from hg19 via the UCSC genome liftover tool (genome.ucsc.deu/cgi-bin/hgLiftOver). For each high-confidence DMR–gene pair, the availability of all CpG sites in the DMR and active expression of the corresponding target gene were verified in the LIHC-TCGA. If no CpGs existed or the gene was not expressed, the replication was said to have failed for that DMR–gene pair. Differential DNA methylation, differential gene expression, and the correlation between DNA methylation and gene expression were evaluated via Wilcox test and Spearman correlation, respectively. A DMR–gene pair was considered to have been replicated in TCGA-LIHC only when there were significant differential methylation, differential gene expression, and consistent (same sign) significant correlation. About the calculation of replication rates of our identified DMR-DEGs in the TCGA-LIHC, it might be unfair to replicate our promoter/enhancer-like DMR-DEGs directly in the TCGA-LIHC cohort, whose DNA methylation levels were profiled by 450k array, especially for DMRs in enhancers, which were rarely covered by 450k array. Thus, DMR-DEGs that failed to replicate were divided into two groups: Group I, replication failure was because no CpG of the DMR was covered in the 450k array; Group II, at least one CpG was available in the 450k array but still failed to replicate the correlated differential methylation and differential expression. Platform-adjusted replicated rates were calculated as: Count_Replicated /(Count_{Type II failure} + Count_Replicated).

In the end, the clinical significance of those high-confidence DMR-DEGs was investigated by way of survival analysis and tumor-stage association analysis. Briefly, the overall survival time (OS) and progression-free survival (PFS) time of TCGA-LIHC samples were retrieved from the integrated TCGA pan-cancer clinical data resource (48). Survival analyses for methylation and gene expression were performed on the basis of the univariate Cox proportional hazards regression model. The Kaplan–Meier method was used to create the survival plots, and the log-rank test was employed to compare the difference between survival curves. The optimal cutoffs of DNA methylation, gene expression, and combination of methylation and expression were determined by minimizing the p-values of log-rank tests. The differences in methylation and gene expression in various tumor stages were compared using ANOVA in R (v.3.5).

For expression detection, 5 × 10⁵ L02 or HepG2 cells were seeded in 12-well plates and allowed to reach 80%–90% confluence. Then, freshly prepared 5-aza-2-deoxycytidine (5-azadC; decitabine) solution was added to the medium at a final concentration of 100 μM. Cells were allowed to grow for 72 h at 37°C with 5% CO2 and then harvested for RNA extraction and qRT-PCR quantification. The cDNA was reverse-transcribed using the iScript™ cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA) according to the manufacturer’s protocol. The real-time PCR was conducted with SYBR Premix Ex Taq (TaKaRa, Kyoto, Japan). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal control for amplification of mRNAs. The comparative C_t method was used to calculate the relative mRNA expression. There were three replicates for each experimental condition (2 cell lines × 2 treatment concentrations).

By WGBS, we obtained the methylation profile of a total of 28,978,826 CpG sites with an average sequencing depth of 12.76 × in 33 pairs of HCC tumor and adjacent non-tumor tissue samples ( Figure 1A , Supplementary Figure S2A and Supplementary Table S1 ). To rescue some CpGs with relatively low depth, smoothed methylation amounts were obtained for all CpGs. We observed a significant methylation difference between tumoral and adjacent tissues and a negative correlation (ρ = -0.49; p = 0.0038) between the average extent of DNA methylation and chronological age in the tumor tissues ( Figure 1B ). This negative correlation also was significant in most genomic regions, including exons, introns, and intergenic regions ( Supplementary Figures S2B–J ). In the PCA plot of all CpGs among all samples, HCCs were clearly separated from paired adjacent tissues, whereas there existed considerable heterogeneities among tumor samples ( Figure 1C ).

Through the combined BWS test, we identified a total of 9,867,700 significant DML between paired HCCs and non-cancerous tissues, including 157,320 hypermethylated DML (hyper-DML) and 9,710,380 hypomethylated DML (hypo-DML). The genomic location annotation of those DML showed that hyper-DML were depleted in the intergenic and intron regions and enriched in other regions, particularly in promoters (defined as -1,500 ~ +500 bp from a TSS) ( Figure 1D ). As for hypo-DML, the pattern was in the opposite direction ( Figure 1D ). Moreover, 47.94% and 43.44% of those hypo-DML were located in the intergenic and intron regions, respectively ( Figure 1D ), which usually were missed in previous arrays or target sequencing-based methylation studies.

After DMR calling from those DML, we identified 608,279 DMRs composed of 6,924 hypermethylated DMRs (hyper-DMRs) and 601,355 hypomethylated DMRs (hypo-DMRs). Promoter annotation of these DMRs revealed 2,882 promoter-like hyper-DMRs and 44,611 promoter-like hypo-DMRs. Of them, 1,569 promoter-like hyper-DMRs and 9,285 promoter-like hypo-DMRs exhibited active promoter-associated histone peaks (H3K4me3) in at least one of those eight liver-related ChIP-seq samples ( Figures 2A, B ) . From their corresponding enhancer activities in tumor and non-tumor liver ChIP-seq samples, we found 3,232 enhancer-like hyper-DMRs and 20,568 enhancer-like hypo-DMRs ( Figures 2C, D ) . Overall, 61.54% of those 6,924 hyper-DMRs were annotated as promoter or enhancer-like regulatory elements ( Supplementary Figure S3A ), while only 4.17% of those 601,355 hypo-DMRs were annotated as promoter or enhancer-like regulatory elements ( Supplementary Figure S3B ). This indicated that hypermethylation events were much less than hypomethylation events but functionally more critical for transcriptional regulation in HCC.

Those promoter-/enhancer-like DMRs were inferred to be “repressed in tumor” ( Figure 2E and Supplementary Figure S3C ) or “activated in tumor” ( Figure 2F ; Supplementary Figure S3D ) according to comparison of their promoter/enhancer activities in tumor and non-tumor samples. In line with the classical negative regulatory relationship between DNA methylation and gene expression, a majority (82.47% and 85.02%, respectively) of those promoter- and enhancer-like hyper-DMRs were recognized as being repressed promoters or enhancers in HCC ( Figures 2G, H ) whereas about half (44.01% and 56.03%, respectively) of those promoter-like or enhancer-like hypo-DMRs were inferred to be repressed promoters or enhancers in HCC ( Figures 2G, H ) , indicating the potential presence of substantial non-classical positive regulation between DNA methylation and histone modification.

The genomic location annotation for all 6,924 hyper-DMRs and 601,355 hypo-DMRs revealed that there were 641 intergenic hyper-DMRs and 250,932 intergenic hypo-DMRs ( Figure S3E ). Quantification of eRNA expression for all intergenic DMRs revealed 36,651 hypo-DMRs and 309 hyper-DMRs with active eRNA expression in our HCC samples. They were recognized as candidate active enhancers in the liver given that eRNA expression is implicated to be a reliable indicator of enhancer activity. The majority of those intergenic active enhancer candidates appeared to be novel. Specifically, 2,378 of them (6.4%) were annotated as active enhancers in at least one of the eight ChIP-seq liver samples, whereas only 6,622 of all 283,631 intergenic DMRs (2.3%) were annotated as active enhancers, an indication of a nearly three-fold enrichment of known enhancers among our intergenic active enhancer candidates. After correlation analysis between DNA methylation and eRNA expression, 4,833 intergenic hypo-DMRs and 23 intergenic hyper-DMRs exhibited a significant negative methylation-eRNA correlation, whereas 2,126 hypo-DMRs and 52 hyper-DMRs displayed non-classical positive methylation-enhancer regulation( Figure S3F ). Only these 7,034 intergenic DMRs with both active eRNA expression and significant methylation–eRNA correlation were defined as intergenic enhancer-like DMRs for downstream analyses.

Aberrant methylation associated genes were determined by integrating DNA methylation with transcriptomic data for our identified promoter-like DMRs, genic enhancer-like DMRs, and intergenic enhancer-like DMRs. Specifically, we found a total of 1,323 potential target genes (i.e., promoter-like DMR-DEGs) for all promoter-like DMRs ( Figure 3A and Supplementary Table S2.1.1 ). As for genic enhancer-like DMRs, we determined 1,751 genes to be their potential targets (i.e., genic enhancer-like DMR-DEGs) ( Figure 3A and Supplementary Table S2.1.3 ). Regarding intergenic enhancer-like DMRs, there were 562 genes (i.e., intergenic enhancer-like DMR-DEGs) that passed the methylation-eRNA–gene triple correlation examination; i.e., these 562 DMR-DEGs displayed simultaneous significant Spearman correlations between DMR methylation and eRNA expression, between eRNA and gene expression, and between DMR methylation and gene expression ( Figure 3A and Supplementary Table S2.1.5 ). Overall, 70.0%, 65.0%, and 79.0% of those identified promoter-like DMR-DEGs, genic enhancer-like DMR-DEGs, and intergenic enhancer-like DMR-DEGs appeared to be negatively regulated by DNA methylation.

The pathway enrichment analysis of negatively correlated DMR-DEGs (i.e., HyperDown and HypoUp) demonstrated deregulated DNA methylations in enhancer-induced activation of genes implicated in the cell cycle, retinoblastoma gene in cancer, DNA replication, and DNA repair ( Figure 3B ) accompanied by the repression of genes in various critical metabolism pathways, including monocarboxylic acid metabolic process, tyrosine metabolism, and carbohydrate metabolic process ( Figure 3C ). Similarly, hypomethylated promoters activated genes implicated in the cell cycle, DNA replication, and DNA repair ( Figure 3D ). On the other hand, 75 genes repressed by promoter hypermethylation failed to be enriched in any pathways, possibly because of the small number of genes, but many of them, such as ST8SIA6-AS1 (49) and GRHL2 (50), were reported to play a suppressor role in multiple cancers including HCC ( Table 2 ). Pathway enrichment analysis of both negatively and positively methylation-correlated DMR-DEGs displayed similar overrepresented pathways ( Supplementary Table S2.2 ).

We identified a set of 611 DMR-DEGs with high confidence through strict screening for genes whose associated DMRs overlapped completely with annotated promoters or enhancers from ChIP-seq ( Table 1 ). Specifically, we discovered 171 high-confidence promoter-like DMR-DEGs ( Supplementary Table S2.1.2 ), 338 high-confidence genic enhancer-like DMR-DEGs ( Supplementary Table S2.1.4 ), and 102 high-confidence intergenic enhancer-like DMR-DEGs ( Supplementary Table S2.1.6 ). Most of the differential DNA methylation in the promoter and genic enhancer regions (70.18% and 73.96%) exhibited a negative correlation with expression of the target gene, but a considerable proportion (41.18%) of those intergenic enhancers showed hypomethylation-associated gene repression.

Literature searching of the identified 56 top differentially expressed high-confidence DMR-DEGs in our HCC sample indicated that 22 of them were implicated in HCC carcinogenesis and the other 15 genes were involved in other types of cancers ( Tables 2 , 3 ). Subsequently, 139/661 high-confidence DMR-DEGs were replicated in the TCGA-LIHC cohort, which consisted of 63 promoter-like DMR-DEGs, 67 genic enhancer-like DMR-DEGs, and 9 intergenic enhancer-like DMR-DEGs ( Table 1 ). Given the limited availability of WGBS-profiled CpGs in the 450k array, particularly the CpGs in non-coding regions, the raw replication rate of DMR-DEGs was modest except for the promoter-like DMR-associated DEGs. Nevertheless, when the platform effect was adjusted, we achieved a considerably higher replication rate (66.32%, 60.36%, and 52.94%) for the above three groups of DMR-DEGs in the TCGA-LIHC cohort ( Table 1 ).

To explore further the clinical significance of our identified 661 potential methylation drivers, we carried out separate survival analyses and tumor stage-association tests based on DNA methylation and transcription data. At the expression level, survival analysis showed that 108 and 82 of those 661 genes were significantly associated with overall survival (OS) and progression-free survival (PFS) ( Figure 4A ), with a high overlap (72 genes) between the two sets of genes. The tumor stage-association test confirmed the significant association between the expression of 140 of the 661 methylation drivers and tumor progression, and 48 of them also were associated with survival ( Figure 4A and Supplementary Table S3.1 ). As for methylation analyses, most of the 661 high confidence DMR-DEGs, except those 139 DMR-DEGs replicated in the TCGA-LIHC cohort, were excluded because of the low coverage in the 450k methylation array. Among those 139 genes, 27, 30, and 18 were highlighted as OS, PFS, and tumor stage-associated methylation biomarkers ( Figure 4B and Supplementary Table S3.2 ). In the same manner, there were extensive overlaps among the three sets of clinically relevant genes ( Figure 4B ). Integration of DNA methylation and gene expression data highlighted six OS-associated, eight PFS-associated, and eight tumor stage-associated DMR-DEGs whose transcription and DNA methylation were both plausible biomarkers ( Figure 4C ). Interestingly, five of those six genes whose expression and methylation were both significant OS biomarkers were more powerful (lower log-rank p-values) prognostic biomarkers when considering both expression and methylation data together ( Figure 4D and Supplementary Table S3.3 ). Similarly, seven of those eight genes whose expression and methylation were both significant PFS biomarkers show higher significance (lower log-rank p-values) when considering both expression and methylation together ( Figure 4E and Supplementary Table S3.3 ). In addition, four genes (CDC20, UCK2, HEATR6, and SLC9A3R2) were concordantly associated with both survival duration and tumor progression ( Figure 4F and Supplementary Tables S3.1, 3.2 ).

For the sake of validation, 15 top differentially expressed STRING protein–protein interaction (PPI) network hub genes of the overrepresented pathways of DMR-DEGs, including cell cycle, DNA repair, and metabolic pathways, and nine genes significantly associated with OS, PFS, and/or tumor stages were selected for in vitro DNA methylation unmasking validation in LO2 and HepG2 cells. After 5-azadC treatment, 20/23 (87.0%) and 15/23 (65.2%) of those selected genes (one gene, named CDC20, belonged to both the pathway hub gene and the clinically relevant gene) showed significant upregulation in LO2 and HepG2, respectively. For instance, four of the five hub genes in the cell cycle and all six hub genes related to the DNA repair pathway showed apparent upregulation in LO2, and a majority of them also presented upregulation in HepG2 after methylation unmasking ( Figures 5A, B ). Likewise, we obtained similar high validation rates for the remaining pathways ( Supplementary Figure S4 ). Furthermore, after 5-azadC treatment in LO2 and HepG2, the qRT-PCR results showed significant upregulation of all seven survival-associated genes and five of the six tumor stage-associated genes ( Figures 5C, D ).