Comprehensive analysis of genomic and immunological profiles in Chinese and Western hepatocellular carcinoma populations

Sample collection

This study includes 100 HCC patients diagnosed from June 2009 to December 2014 in the West China Hospital, Sichuan University. The tissue samples in this study include tumor tissue and matched liver tissue adjacent to the tumor. The application and use of samples were reviewed and approved by the Institutional Ethics Committee of the West China Hospital. In addition, this study include 175 patients from the TCGA (http://cancergenome.nih.gov/) database. HBV-related HCC from the West China Hospital was abbreviated as WCH-HBV-HCC; HBV-related HCC in the Western population of TCGA database was abbreviated as TCGA-HBV-HCC; NonHBV-related HCC in the West China Hospital was abbreviated as WCH-NonHBV-HCC; NonHBV-related HCC in the Western population from the TCGA database was only included alcoholic HCC, abbreviated as TCGA-Alcol-HCC. All HBV-related HCC patients were positive for hepatitis B surface antigen (HBsAg). The main causes of 19 nonHBV-related HCC patients in the West China Hospital are as follows: 3 metabolic syndrome related diseases (diabetes and obesity), 2 alcoholic, and the remaining 14 patients have no obvious cirrhosis and no clear cause. Tumor staging is based on the eighth edition of the American Cancer Society (AJCC)-TNM staging. The degree of tumor differentiation according to Edmondson-Steiner standard is divided into I-IV grade. I-II grade is defined as high differentiation and III-IV grade is defined as low differentiation.

Immunohistochemical staining

The tissue microarray is constructed by the Xinchao Biotechnology (Shanghai, China). The steps of the immunohistochemical (IHC) staining were described previously [20]. Details of the primary antibodies used in IHC staining is shown in the Supplementary Table 1.

In the WCH-HBV-HCC and WCH-NonHBV-HCC groups, the number of three types of immune cells were evaluated by IHC staining including pan-T cells (CD3), B cells (CD20) and cytotoxic T cells (CD8). Before counting immune cells, five representative areas were selected (each area is 1mm²) at 100 × magnification, and then the number of positive immune cells was counted at 200 × magnification. The average of the five areas is the number of positive immune cells. The values of other immune markers were recorded as cells/mm². Methods for analysis of tumor-infiltrating immune cells for samples from the TCGA-HBV-HCC and TCGA-Alcol-HCC groups was shown in the Supplementary Materials and Methods.

Whole exome sequencing

After DNA extraction (UPure Tissue DNA Kit, BioBase, Sichuan, CN), we performed DNA quality control in three aspects: (1) Evaluation of DNA degradation by agarose gel electrophoresis; (2) Agilent Bioanalyzer 2100 (Agilent, USA) detects the size of DNA sample fragments; (3) Thermo Scientific™ NanoDrop™ (Thermo Fisher Scientific, USA) tests the purity of DNA samples (OD260/280 and OD230/260).

The DNA was fragmented using the ultrasonic interrupter Bioruptor® (Diagenode, Belgium). After connecting the sequencing adaptor, the Illumina sequencing library was constructed utilizing Nextera DNA Library Prep Kit (Illumina, USA). For exon capture, Agilent Sureselect™ Human All Exon V6 kit (Agilent, USA) was used for liquid phase hybrid capture of the exon portion of DNA. After building the high-throughput sequencing library, the sequencing was performed on the platform Illumina Novaseq™ 6000 (Illumina, USA) with the sequencing mode of PE150. After removing the interference from PCR duplicates, the average sequencing depth of cancer tissues was 200 ×, and the average sequencing depth of adjacent tissues was 100 ×.

Genomic analysis

We firstly used the FASTQ software to remove the joint sequences and low-quality sequences in the original data to obtain the "clean data" in the format of FASTQ. The BWA software was used to perform reference genome alignment (hg19 version) with the reads contained in paired FASTQ files, and the Picard software was utilized to identify and mark duplicate reads from BAM file. Next, GATK4.0 software was used to identify the mutations of SNP and Indel. In this study, when copy number variation (CNV) was detected, the sequencing depth was firstly corrected based on the GC content of the sequence, and finally the presence of CNV was determined after normalization of the data by clustering analysis. Tools used for CNV analysis includes GISTIC 2.0 [21] and CNVkit [22]. The significance of gene mutation was identified by MutsigCV [23] and visualized by maftools and ggplot2. The non-negative matrix factorization (NMF) method was based on the previous study [24]. The identified signatures were compared with the COSMIC signatures.

Analysis of data from TCGA

The TCGA genomic and transcriptome sequencing data was obtained from the Genomic Data Commons (GDC) Data Portal (https://portal.gdc.cancer.gov/). We utilized the same analytic methods of genomic analysis for patients in the WCH and TCGA cohorts. The analytic methods of tumor-infiltrating immune cells in the TCGA-HBV-HCC and TCGA-Alcol-HCC cohorts were shown in the Supplementary Materials and Methods in detail.

Statistical analysis

Categorical variables were expressed as number (%), and the statistical test method was Chi-square test or Fisher^, exact test. Continuous variables were expressed as mean ± standard deviation, and the statistical method was T-test or Kruskal-Wallis test. The FDR-corrected q values were used to compare the differences in frequency of the common gene mutation or CNV number between two groups, and the calculation method was the Benjamini-Hochberg method. Data analysis was carried out by R software 3.4.3.

Sample and clinical information

The clinicopathologic information of HCC patients from the West China Hospital is shown in Supplementary Table 2. Compared to patients in the WCH-HBV-HCC group, patients in the WCH-NonHBV-HCC group were older (P = 0.002), while the other clinicopathologic features were not significantly different between the two groups (all P values > 0.05). For patients in the TCGA group, patients with alcoholic HCC (TCGA-Alcol-HCC group) were older than those with HBV-related HCC (TCGA-HBV-HCC group) (P <0.001). In addition, compared to patients from the TCGA-HBV-HCC group, cases in the TCGA-Alcol-HCC group had higher proportion of white people (P <0.001), later TNM stage (P <0.001), better tumor differentiation (P <0.001), and lower proportion of microvascular and macrovascular invasion (P = 0.032). The remaining characteristics had no significant difference between the two groups (all P values > 0.05) (Supplementary Table 3).

Somatic copy number variation

In this study, GISTIC 2.0 was used to identify copy number changes for patients in the HBV and non-HBV groups, including copy number amplification and deletion. We found that, there was a significant difference in total CNV levels for patients in the WCH-HBV-HCC group (Supplementary Figure 1B) and WCH-nonHBV-HCC group (Supplementary Figure 1C).

For patients in the WCH-HBV-HCC group, copy number amplification (CNV-Amp) was found mainly on chromosomes 1 and 8, and other locations including 4p, 4q, 5p, 6p, 7p, 7q, 11q, 12q, 13q, 14q, 15q, 17p, 17q, 19p and 19q, etc. Copy number deletion (CNV-Del) was observed predominantly on chromosomes 7q, 8p, 15q, and other locations including 1p, 1q, 2q, 3q, 4p, 4q, 5q, 6q, 9p, 9q, 10q, 11p, 11q, 16p and 17q, etc. Among WCH-NonHBV-HCC patients, CNV-Amp was found to mainly locate at 1q, 8q, 10q, 11p, and 19p; CNV-Del was found to mainly locate at 8p, 15q, and 9p (Figure 1). In summary, the CNV levels in the WCH-HBV-HCC group and WCH-NonHBV-HCC group are significantly distinct.

Firstly, we compared the CNV levels of patients in the WCH-HBV-HCC group and the TCGA-HBV-HCC group. We found that in the above two groups, the locations and frequencies of CNV-Amp and CNV-Del were obviously different between two groups (Supplementary Figure 2A, 2B). For example, cases in the WCH-HBV-HCC group had CNV-Amp in 4p, 4q, 5p, 6p and 14q and CNV-Del in 14q, 15q, 16p and 21q, while CNV in these locations were not significantly observed in the TCGA-HBV-HCC group. Besides, significant difference was also observed between the WCH-NonHBV-HCC group and the TCGA-Alcol-HCC group (Supplementary Figure 2C, 2D).

Additionally, patients in the WCH-HBV-HCC group and the TCGA-HBV-HCC group were compared at the gene level (genes affected by CNV). The number of unique CNV-Amp genes in the WCH-HBV-HCC group and TCGA-HBV-HCC group were 668 and 771, respectively, while the number of CNV-Amp genes shared by the two was 184. The number of unique CNV-Del genes of the above two groups is 2487 and 236, respectively, while the number of CNV-Del genes shared by the two groups is 684. Supplementary Figure 3 shows parts of the CNV-affected genes that are common or unique in the WCH-HBV-HCC and TCGA-HBV-HCC groups (Supplementary Figure 3A, 3B), or WCH-NonHBV-HCC and TCGA-Alcol-HCC groups (Supplementary Figure 3C, 3D). The shared CNV-Amp genes include CYC1, ATP4B, SCCPDH and KISS1, etc. The shared CNV-Del genes include ACADVL and FAM138A, etc. In addition to the above-mentioned shared CNV genes, there are a large number of CNV-affected genes that only observed in their respective cohorts. For example, CNV-Amp of MUC16, ZDHHC1, DUX4, NBPF20, PPIAL4C, REXO1L2P, RASA4 and CCND1 was only found in the WCH-HBV-HCC group, and CNV-Amp of ADCY8, TPCN2, FOXK2 and CAPZA2 was only shown in the TCGA-HBV-HCC group. In addition, in line with expectations, due to differences in pathogenic causes, the locations of CNV (Supplementary Figure 2C, 2D) and the genes affected by CNV (Supplementary Figure 3C, 3D) for patients in the WCH-NonHBV-HCC and TCGA-Alcol-HCC cohorts are significantly heterogeneous.

In addition to the genes shown in Supplementary Figure 3, there are a large number of common genes affected by CNV that appear separately in the WCH-HBV-HCC group and the TCGA-HBV-HCC group (Supplementary Figures 1–4). For example, CNV-Amp of genes such as FGF4, FGF19, ETV1, PMS2, RPTOR, CARD11, ERCC5 and RAC1 appeared only in the WCH-HBV-HCC group, and CNV-Amp of genes such as MET, MYC, AKT3, AGO2, BRD4 and IKBKE appeared only in the TCGA-HBV-HCC group. In addition, we found that CNV-Del of genes such as CDKN2A, AXIN1, ACVR1, PTPRD and ATM only appeared in the WCH-HBV-HCC group, while CNV-Del of genes such as FGFR2, STK11, TCF3, and KMT5A were only shown in the TCGA-HBV-HCC group.

Somatic mutation and tumor mutation burden

There are 4.5 (median) somatic mutations detected per million bases in WCH-HBV-HCC group, and the number is 5 (median) in WCH-NonHBV-HCC group. As shown in Figure 2, the most common mutational genes in patients from the WCH-HBV-HCC group are TP53, TTN, MUC4, CTNNB1, CDC27 and MUC16, etc. The most common mutations in patients from the WCH-NonHBV-HCC group are MUC16, TP53, TTN, ASTN1, ARID1A and HMCN1, etc. The most common mutant genes in cases from the TCGA-HBV-HCC group are TP53, CTNNB1, TTN, MUC16, AXIN1 and MUC4, etc. The most commonly mutated genes for patients in the TCGA-Alcol-HCC group are TP53, CTNNB1, TTN, PCLO, ALB and MUC16, etc. We found that more than half (51%) of patients in WCH-HBV-HCC cohort had TP53 mutations, while CTNNB1 mutations accounted for only 28% (ranking fourth). Other highly-mutated bur rarely-reported genes include TTN, MUC4, CDC27, MUC16, PCLO and UBA52, etc. The most frequently mutated genes in WCH-NonHBV-HCC patients are MUC16 and TP53, and only 4 patients in this group had CTNNB1 mutation. In contrast, in the TCGA-HBV-HCC group and the TCGA-Alcol-HCC groups, TP53 and CTNNB1 mutations are the most common genetic mutations. The sum of patients with TP53 or CTNNB1 mutation exceeds half of the individual cohort. Most of the samples do not have TP53 and CTNNB1 mutations simultaneously, which indicates that these two mutations represent two distinct types of tumors.

Similar to CNV, the main driver mutations showed significant heterogeneity in the four cohorts including WCH-HBV-HCC, WCH-NonHBV-HCC, TCGA-HBV-HCC and TCGA-Alcol-HCC groups (Figure 3 and Supplementary Figure 4). Gene with q value <0.1 was defined as significantly mutated gene. Figure 3E shows the significantly mutated genes appearing in the above groups. The genes presented in Figure 3E are those with q values <0.1 in at least one group. This study compared the significantly mutated genes in the above four groups and found that TP53 and UBA52 are the significantly mutated genes shared in the WCH-HBV-HCC group and WCH-NonHBV-HCC group (Figure 3F); TP53, CTNNB1 and ALB are significantly mutated genes shared in the TCGA-HBV-HCC group and the TCGA-Alcol-HCC group (Figure 3G); TP53, CTNNB1 and AXIN1 are the significantly mutated genes shared in the WCH-HBV-HCC group and the TCGA-HBV-HCC group (Figure 3H); TP53 is a significantly mutated gene shared in the WCH-NonHBV-HCC group and the TCGA-Alcol-HCC group (Figure 3I). Consequently, TP53 is a significantly mutated gene shared by the four groups of patients in this study, and the role of UBA52 mutations in HCC has not been studied (Figure 3F). In addition to these shared genes, a large number of unique mutations were observed in each groups (Supplementary Figure 4 and Figure 3E–3I).

Mutation signature and mutation spectrum analysis

There are 6 types of point mutations: C>A/G>T, C>G/G>C, C>T/G>A, T>A/A>T, T>C/A>G, T>G/A>C. We analyzed the mutation spectrum of patients in the WCH-HBV-HCC group and the WCH-NonHBV-HCC group (Supplementary Figure 5), and found that the two most commonly mutation types were C>T and T>A in the WCH-HBV-HCC group, followed by T>C and C>A mutations. In contrast, in the WCH-NonHBV-HCC group, the frequencies of T>C and C>A mutations were significantly lower than the C>T and T>A mutations.

Based on results (cophenetic curves) of the Nonnegative Matrix Factorization (NMF) method (Supplementary Figures 6, 7), we found that the optimal number of clusters in both groups is 3. Finally, mutation signature analysis of 96 substitution patterns identified 3 signatures in the WCH-HBV-HCC group and WCH-NonHBV-HCC group, respectively (Figure 2A, 2B). Signature A of the WCH-HBV-HCC group and WCH-NonHBV-HCC group is a mutation signature that has been described and verified in previous studies, namely Signature 22 (characterized by dominant T>A mutations) in the COSMIC database, which is a characteristic mutation of HCC related to aristolochic acid. Signature B and Signature C in the WCH-HBV-HCC group are newly discovered mutation signatures (Figure 2A). Signature B (C>T) of the WCH-HBV-HCC group shows some similarity to the previously described Signature 1 (characterized by dominant C>T mutations) and 6 (C>T mutation; related to the loss of DNA mismatch repair function), 14 (C>T and C>A), 15 (C>T), 19 (C>T) and 23 (C>T). Signature C of WCH-HBV-HCC group has part of characteristics of COSMIC Signature 5 (T>C), 16 (T>C), 3 (shows all 6 types of mutation), 12 (T>C) and other signatures (Supplementary Figure 8A). Signature B and Signature C of the WCH-NonHBV-HCC group are also newly discovered mutation signatures (Figure 2B). Signature C (C>T) of the WCH-NonHBV-HCC group shows some mutation features of COSMIC Signature 1 (C>T) and 6 (C>T). Signature B of the WCH-NonHBV-HCC group has some similarity to COSMIC Signature 5 (T>C), 16 (T>C), 3 (all 6 mutation types appear) 4 (C>A; related to smoking) and 8 (C>A) (Supplementary Figure 8B). In addition, previous literatures reported that COSMIC Signature 24 (C>A) is a predominant signature of aflatoxin-related HCC. In this study, Signature B of WCH-HBV-HCC group and Signature C of WCH-NonHBV-HCC group have some mutation characteristics of Signature 24.

Figure 2C–2E, 2F–2H display the distribution of mutation signatures (Signature A, B and C) of all samples in the WCH-HBV-HCC group and WCH-NonHBV-HCC group, respectively. We found that samples in both groups showed obvious heterogeneity, and the dominant signature was distinct in different samples. Meanwhile, we clustered the HCC samples of the WCH-HBV-HCC and WCH-NonHBV-HCC groups according to the 30 mutation signatures identified in the COSMIC database (Supplementary Figures 9, 10). We observed that, except for Signature A (similar to COSMIC Signature 22), there were no signatures in the WCH group showed a strong similarity to the previously described signatures.

Immune cell infiltration in different types of HCC

We analyzed the immune cell infiltration in the TCGA-HBV-HCC group and the TCGA-Alcol-HCC group based on the mRNA sequencing data (Supplementary Figure 11). Based on the specifically expressed gene set of immune cells related to innate and adaptive immune responses, using GSVA method, we calculated the relative expression of immune cells, angiogenesis and antigen-presenting machinery (APM) in the above two groups. The results showed that cases in the TCGA-HBV-HCC group and the TCGA-Alcol-HCC group had numerous T cell infiltration. Compared to cases in the TCGA-Alcol-HCC group, patients in the TCGA-HBV-HCC group had more infiltration (without statistically significance) of helper T cells (96.12% vs. 87.30%), overall T cells (68.93% vs. 57.14%) and B cells (9.71% vs. 1.59%). In contrast, the distribution of the remaining immune cells between the two groups is more similar. The above results indicate that, compared to gene alterations, the immune cell infiltration in the TCGA-HBV-HCC group and the TCGA-Alcol-HCC group is less heterogeneous. It is worth noting that a high frequency of samples in both groups have infiltration of CD8+ T cells and other types of T cells, which provides the possibility of immunotherapy.

Association of gene mutation with immune cell infiltration in the TCGA cohorts

Using the GSVA algorithm, we calculated the distribution of 24 types of immune cells for patients in the TCGA-HBV-HCC group (Figure 4A) and TCGA-Alcol-HCC group (Figure 4B). Meanwhile, after hierarchical clustering analysis based on the GSVA score, we divided patients in the above two groups into three subgroups with distinct immune cell infiltration status (Figure 4). The immune-HIGH subgroup has the highest GSVA score and more immune cell infiltration and the GSVA score and number of tumor-infiltrating immune cell in the immune-LOW subgroup is the lowest. The GSVA score and immune cell number in the immune-MIX subgroup is between the immune-HIGH and immune-LOW groups. As shown in Figure 4, there is no significant relationship between tumor mutation burden and immune cell infiltration for patients in the TCGA-HBV-HCC and TCGA-Alcol-HCC groups. Besides, we analyzed the relationship between gene mutations and immune cell infiltration. We found that, compared to the total number of patients in the immune-HIGH and immune-MIX group, more cases in the immune-LOW subgroup had TP53 mutation (P = 0.007) (Figure 4A). In contrast, in the TCGA-Alcol-HCC group, compared to all patients in the Immune-HIGH and Immune-MIX subgroups, cases in the Immune-LOW subgroup had similar frequency of mutations in TP53, CTNNB1 and other genes (all P values> 0.05) (Figure 4B).

Association of gene mutation with immune cell infiltration in the WCH-HBV-HCC group

For patients in the WCH group, all tumor samples were performed immunohistochemical staining. Immune cell markers include CD3, CD8 and CD20. A representative image of the above IHC staining is shown in Supplementary Figure 12. The hematoxylin and eosin (HE) staining shows that cases in the WCH-NonHBV-HCC had distinct histopathological features (Supplementary Figure 13). Owing to the limited sample number and heterogeneous composition in the WCH-NonHBV-HCC group, in this section, we only analyzed the association of gene mutation and tumor-infiltrating immune cell in the WCH-HBV-HCC group (Figure 5). We divided the WCH-HBV-HCC group into Immune-HIGH, Immune-MIX, and Immune-LOW subgroups (Figure 5A). Similar to the TCGA cohorts, the number of most types of immune cells including CD8+ T cells, CD3+ T cells and CD20+ B cells is the highest in the immune-HIGH subgroup. Consistent to patients in the TCGA-HBV-HCC group, the tumor mutation burden in the WCH-HBV-HCC group was not significantly associated with the infiltration of immune cells in the tumor. Among patients in the WCH-HBV-HCC group, we verified the findings in the TCGA-HBV-HCC group. Finally, we found that the Immune-LOW subgroup had more TP53 mutation than cases in the immune-HIGH and immune-MIX subgroups (P = 0.049).

This study also analyzed the relationship between TP53 mutation types and immune grouping (Figure 5B). The main mutation types of TP53 in the WCH-HBV-HCC group include Missense mutation, Nonsense mutation, Frame shift deletion and Splice mutation. Figure 5A clustered the immune cells based on the above mutation types and find that, for patients in the WCH-HBV-HCC group, all frame-shift mutations, nonsense mutations, and parts of missense mutations of TP53 were associated with the decrease in the number of tumor-infiltrating immune cells.

Association of CNV with immune cell infiltration

We explored the relationship between the immune infiltration score (ESTIMATE immune score) and gene copy number alterations. We observed that higher number of copy number amplifications correlated with lower ESTIMATE immune infiltrate scores in both TCGA-HBV-HCC (Figure 6A) and TCGA-Alcol-HCC (Figure 6C) groups. In addition, higher number of copy number losses were also associated with lower ESTIMATE immune scores in both TCGA-HBV-HCC (Figure 6B) and TCGA-Alcol-HCC (Figure 6D) groups (P < 0.05, Pearson’s correlation). To identify genes with CNVs that were significantly associated with immune cell infiltrate, we evaluated the relationship between the most commonly CNV-affected genes of HCC (VEGFA, FGF3/4/19, CCND1, AXIN1, CDKN2A, CDKN2B, IRF2, MAP2K3, PTEN and RB1) and ESTIMATE immune infiltrate scores. However, none of the above genes^, CNV were related to the immune infiltration scores. In melanoma, PTEN loss was found to be associated with the immune cell infiltration status, while on such associations were observed in HCC (TCGA-HBV-HCC, WCH-HBV-HCC and TCGA-Alcol-HCC groups) (Figure 6E).