The transcriptomic (n = 39) and methylation (n = 70) profile data, as well as the clinical and pathologic characteristics of Peruvian patients diagnosed with HCC and programmed for liver resection, within the time frame of 2006 and 2016 at the INEN, were obtained from a previous study and downloaded from the Gene Expression Omnibus (GEO) database.15 In brief, none of the samples collected came from patients that had undergone chemotherapy before liver resection, in accordance with the recommendations of the European Association for the Study of the Liver (EASL).17 The samples were collected in pairs; an HCC tissue sample and an adjacent nontumor liver tissue (NTL) sample extracted from the surgical margin (≥2 cm). Approximately 50 mg of each tissue type (HCC and NTL) were flash frozen and stored at −150 °C in the Tumor Tissue Bank of the INEN before analysis. The variables considered in the present study included demographic, laboratory, histologic, and HBV infection data. Table S1 shows the details of each variable.

The patients were classified according to the survival prognosis nomogram for tumors >10 cm.18 The preoperative variables of that nomogram are: i) presence of cirrhosis, ii) multiple tumor lesions, iii) macroscopic vascular invasion, and iv) spontaneous tumor rupture. From the final score, the samples were classified into three survival groups: G1, with a high probability of long-term survival (0 points); G2, with a 60% probability of recurrence within 12 months (between one and 100 points); and G3, with the highest postoperative mortality (>100 points).

HBV DNA detection in tissues was carried out as previously described.19 In short, genomic DNA was extracted using the phenol-chloroform method, and quantification was done, using the Qubit® (Invitrogen) dsDNA kit. All samples were examined for detecting the presence of HBV DNA, utilizing the PCR QX200™ Droplet Digital™ (Bio-Rad) system. HBV sequencing was performed using the BigDye™ Terminator v3.1 (Applied Biosystems) sequencing kit and HBV phylogeny was determined using the two-parameter Kimura model and the neighbor joining method in version 7 of the Molecular Evolutionary Genetics Analysis software (MEGA, University Park, PENN, USA) The sequencing data can be downloaded from the European Bioinformatics Institute (EBI) site, with the access code: PRJEB21100 (https://www.ebi.ac.uk/ena/data/view/PRJEB21100).

The transcriptome profile was analyzed from a group of Peruvian HCC samples obtained from the Gene Expression Omnibus database (GEO, Bethesda, MD, USA), with codes GSE136247 and GSE111580.15 The transcriptome of those samples was obtained using GeneChip™ Human Transcriptome Array 2.0 microarrays and the GeneChip™ WT PLUS Reagent Kit (Applied Biosystems, Waltham, MA, USA). The transcriptomic data were normalized through the Robust Multichip Average (RMA) and the lot effect was corrected using the ComBat algorithm of the limma package in R (R language).20 The matrix (genes x samples) was then collapsed using the collapse_genes-09 (Fred’s Softwares, Toulouse, France) algorithm. The differentially expressed genes (DEGs) were obtained, using the R limma package, with a corrected p value (q value) of 0.05, following the recommended workflow.20 In addition, the 5000 most variable genes in the HCC tissue were identified, using standard deviation. Heat maps based on unsupervised hierarchical clustering were created, using the pheatmap R package. Identification of the metabolic pathways of interest was carried out using the Sample Enrichment Score (SES) algorithm and the KEGG, Reactome (reactome_2018803), Gene Ontology-Cellular Component (C5_CC), and MSigDB Hallmarks (H) databases, compiled by Fred’s Softwares.21 The source codes and databases of Fred’s Software are available at: https://sites.google.com/site/fredsoftwares/products.

The methylation profile of the Peruvian HCC samples, with codes GSE136319 and GSE136380, were obtained from the GEO database. The methylome of the samples was obtained, using the Infinium HumanMethylation450 K and MethylationEPIC BeadChip (Illumina, San Diego, CA, USA) microarrays. The differentially methylated positions (DMPs) were identified, using the minfi R package, with a q value of 0.05.15 The criteria considered for the methylation analysis were: mean Δβ between HCCs and NTL >5% (for the hypermethylated positions) or ≤ 5% (for the hypomethylated positions). The 5000 positions with the most variable methylation were identified through calculating the standard deviation and the classification, regarding the association with “CpG islands” and “gene structure”, was carried out according to the Illumina definition. The heat maps based on unsupervised hierarchical clustering were created using the pheatmap R package. The SES algorithm was utilized to identify the metabolic pathways of interest and the genes associated with the CpGs of interest, according to the Illumina definition, were used.

The multiple factor analyses (MFAs) and principal component analyses (PCAs) were carried out in R, using the FactoMineR package.22 Variables with incomplete data were handled with the imputeMFA function of the missMDA R package. In brief, as recommended in the FactoMineR user guide, an MFA was applied, given that our data had numerical and categorical variables. The MFA was done, according to the recommendations established in the package. All the groups were considered relevant variables and no supplementary data were incorporated in the analysis. In addition, the t-test and chi-square test were done in R. Statistical significance was unilaterally or bilaterally evaluated, depending on the circumstances. All the figures presented herein were created in R and GraphPad Prism v8 (GraphPad Software, Boston, MA, USA).

All the samples included in in the present study came from volunteer patients that accepted the use of their samples and data (in anonymity) for research. The project (code INEN 10-05) was reviewed and approved by the Institutional Committee on Ethics in Research of the INEN.

Seventy-four Peruvian patients with HCC that underwent liver resection at the INEN were included in the study. Mean patient age was 45.9 years. In accordance with the median age (44 years), the patients ≤44 years of age were classified as adolescents and young adults (AYAs) and patients >44 years of age were classified as middle or older age individuals (MOAs). The AYA group accounted for more than 50% (n = 39) of the study population with HCC, with a mean age of 27 years, and the MOA group (47.3%; n = 35) had a mean age of 67 years. The male-to-female ratio for the entire patient cohort was two males for each female (2:1) but was double (4:1) in the MOA group. Mean tumor diameter was 14.4 cm (AYAs 14.2 cm and MOAs 14.6 cm; p > 0.05) and the large majority of patients (∼80%) presented with a tumor > 10 cm in diameter. The survival prognosis nomogram for tumors >10 cm,18 was retrospectively applied to 36 patients (21 AYAs and 15 MOAs), and 83% of those patients were in the G1 (n = 17) and G2 (n = 13) groups, both of which had a higher survival probability, compared with the G3 (n = 6) group, in which all the patients died before reaching the three-year postoperative period (Fig. 1a). The Hepatitis B surface antigen (HBsAg) was detected in 53% of the patients (72% of AYAs and 31% of MOAs), whereas HBV DNA was detected in 86% of the study population with HCC; the only sub-genotype identified was F1b. Lastly, the AYA and MOA groups accounted for 92 and 80% (p > 0.05) of the cases with viral DNA, respectively (Table 1).

Figure 1.

Evaluation of the clinical and demographic variables in Peruvian patients with HCC. (a) G1 (n = 17), G2 (n = 13), and G3 (n = 6) survival curves of 36 patients, classified according to nomogram. As a result of the MFA, the main components are shown in the graph: (b) the contribution of each group of variables, and (c) the most relevant variables within each group, for the 2 most relevant dimensions (Dim. 1 and 2). Lab (blue); Demog (green and red); HBV (light blue); Histo (yellow and grey); Tumor (brown and fuchsia). Qualitative (ql) and quantitative (qt) data. The PCA of (d) gene expression and (e) DNA methylation data of Peruvian patients, according to HBV infection, age (AYAs and MOAs), and tissue type (HCC and NTL).

(0.64MB).

For the MFA, the variables were organized into eight groups: clinical laboratory data (Lab), qualitative demographics (Demog ql), quantitative demographics (Demog qt), cases associated with HBV infection, qualitative histologic (Histo ql) variables, quantitative histologic (Histo qt) variables, qualitative tumor-associated (Tumor ql) variables, and quantitative tumor-associated (Tumor qt) variables (supplementary Table S1). The MFA, utilizing all the clinical and demographic data of the patients with HCC, showed that the groups of the Demog ql, Histo ql, Histo qt, Lab, and Tumor qt variables were the greater contributing parameters (Fig. 1b). The variables within those groups with higher relevance were age, survival, number of mitotic cells (mitosis), AFP, tumor diameter, total bilirubin, and albumin, according to the MFA (Fig. 1c).

Due to the high HBV infection rate, a PCA was done to determine the effect of HBV on gene expression and methylation. In the two settings, there was a higher number of AYA patients infected with HBV, showing that HBV acts as an age-related confounding factor (Fig. 1d, e).

The transcriptomic profiles in HCC and NTL from a total of 39 patients (19 AYAs and 20 MOAs) were obtained from previously published data.15 The evaluation of the 5000 most variable genes in HCC showed no significant difference between the two age groups (p > 0.05) (Fig. 2a). However, an unsupervised hierarchical cluster analysis of the 5000 most variable genes in the HCC tissue showed the genes were divided into three large groups (C1, C2, and C3). The genes in C1 were mainly associated with expression regulation through histone modification, whereas the majority of genes in C2 participated in the complement cascade and coagulation cascade, and the genes in C3 participated in the immune response. Other less representative pathways were hormone regulation (C1), retinol and lipid metabolism (C1), and DNA repair (C2) (Fig. 2b).

Figure 2.

Analysis of the transcriptomic and DNA methylation profiles of Peruvian patients with HCC. 5000 genes (a,b) and 5000 most variable CpGs (c,d) were identified. Utilizing (a) the expression (AYA = 19, MOA = 20) and (c) DNA methylation levels (AYA = 39, MOA = 35) heat maps were created based on an unsupervised hierarchical cluster. In (a), three large gene groups (C1, C2, and C3) stand out. The pathways with greater significance and fold enrichment (b) of the 5000 genes, and (d) of the 5000 most variable CpGs are represented in a dispersion graph. Fold enrichment is the number of genes from the list of interest that are found on the pathway evaluated.

(1.35MB).

The methylation profiles of HCC and NTL from 70 patients (36 AYAs and 34 MOAs) were obtained from a previous study.15 Similar to the transcriptome analysis, the 5000 most variable CpGs were identified and unsupervised hierarchical clustering was performed. In contrast to the NTL, there was an overall increase in DNA methylation (hypermethylation) in the Peruvian HCC, with greater intensity in the younger patients (Fig. 2c). Those CpGs were also mainly associated with the genes involved in cell differentiation, with methylation marks of lysine residue at N-terminal position 27 of the H3 histone protein (H3K27 [me3]) and targets of the Polycomb Repressive Complex 2 (PRC2) (Fig. 2d).

To better understand CpG DNA methylation, a heat map was created on the positively and negatively regulated genes.15 The regions of gene structure and CpG islands were considered independently. The CpGs associated with negatively regulated genes (the majority hypermethylated) were mainly detected on the CpG islands, within the TSS 200 and TSS 1500 promotor regions, whereas the CpGs associated with positively regulated genes (the majority hypomethylated) were located at regions unrelated to the CpG islands in the gene body (Fig. 3a, b). Genomic integration in HCC also revealed a predominance of the pathways related to cell metabolism, but pathways related to cell differentiation and the immune system stood out, as well.

Figure 3.

Integrated analysis of the transcriptome and DNA methylation profiles of Peruvian patients with HCC. Heat maps based on an unsupervised hierarchal cluster to show the DNA methylation levels (β value) of (a) the negatively and (b) positively regulated differentially expressed genes (DEGs) in the AYA (n = 39, upper panel) and MOA (n = 35, lower panel) groups. (a,b) The left dendrogram shows the annotation of the CpG island and gene structure. (c) The fold enrichment pathways of greater significance associated with the integrated analysis of gene expression and DNA methylation are illustrated, utilizing a dispersion graph. The fold enrichment is the number of genes on the list of interest that are found on the pathway evaluated.

(1.57MB).

Lastly, the differences between the groups defined according to the survival prognosis nomogram for tumors >10 cm, at the transcriptomic and DNA methylation levels, were explored. Only 36 patients (17 for the transcriptomic profile and 36 for the DNA methylation profile) had the preoperative variables necessary for being evaluated. A total of 17 patients obtained 0 points, and so were classified as G1, whereas 13 patients obtained between 1–100 points and 6 obtained >100 points, and so were classified as G2 and G3, respectively (Table 1). According to that classification, 47% of the patients had a high probability of long-term survival, greater in those ≤44 years of age (12 AYAs and 5 MOAs), 36% had a high probability of recurrence within the first 12 months (7 AYAs and 6 MOAs), and 17% had a high probability of postoperative death (2 AYAs and 4 MOAs). Despite those clinical differences, the analyses showed no significant differences at the transcriptomic or DNA methylation levels, between the three groups; thus, those analyses were not studied in more detail. Nevertheless, more than 50% of the patients ≤ 44 years of age (57%) had a high probability of long-term survival.