A Novel Immune Classification Reveals Distinct Immune Escape Mechanism and Genomic Alterations: Implications for Immunotherapy in Hepatocellular Carcinoma

doi:10.21203/rs.3.rs-100244/v1

Download PDF

Research

A Novel Immune Classification Reveals Distinct Immune Escape Mechanism and Genomic Alterations: Implications for Immunotherapy in Hepatocellular Carcinoma

https://doi.org/10.21203/rs.3.rs-100244/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 06 Jan, 2021

Read the published version in Journal of Translational Medicine →

You are reading this latest preprint version

Background: The tumor immunological microenvironment (TIME) has a prominent impact on prognosis and immunotherapy. However, the heterogeneous TIME and the mechanisms by which TIME affects immunotherapy have not been elucidated in hepatocellular carcinoma (HCC).

Methods: A total of 2195 eligible HCC patients from TCGA and GEO database were collected. We comprehensively explored the different heterogeneous TIME phenotypes and its clinical significance. The potential immune escape mechanisms and what genomic alterations may drive the formation of different phenotypes were further investigated.

Results: We identified three phenotypes in HCC: TIME-1, the “immune-deficiency” phenotype, with immune cell depletion and proliferation; TIME-2, the “immune-suppressed” phenotype, with being in immunosuppressive states; TIME-3, the “immune-activated phenotype”, with abundant leukocytes infiltration and immune activation. The prognosis and sensitivity to both sorafenib and immunotherapy differed among the three phenotypes. We also underlined the potential immune escape mechanisms: lack of leukocytes and defective tumor antigen presentation capacity in TIME-1, increased immunosuppressive cells in TIME-2, and rich in immunoinhibitory molecules in TIME-3. The different phenotypes also demonstrated specific genomic events: TIME-1 characterized by TP53, CDKN2A, CTNNB1, AXIN1 and FOXD4 alterations; TIME-2 characterized by significant alteration patterns in the PI3K pathway; TIME-3 characterized by ARID1A mutation. Besides, the TIME index (TI) was proposed to quantify TIME infiltration pattern, and it was a superior prognostic and immunotherapy predictor. A pipeline was developed to classify single patient into one of these three subtypes and calculated the TI.

Conclusions: We identified three TIME phenotypes with different clinical outcomes, immune escape mechanisms and genomic alterations in HCC, which could present strategies for improving the efficacy of immunotherapy. TI as a novel prognostic and immunotherapeutic signature that could guide personalized immunotherapy and clinical management of HCC.

Translational Medicine

Hepatocellular carcinoma

Immunotherapy

Tumor immunological microenvironment

Immune escape

Molecular subtype

Hepatocellular carcinoma (HCC) is the dominant histologic type of primary liver cancer, with a high incidence and mortality rate [1]. Although there are various therapeutic modalities for HCC, including surgical resection, chemotherapy, radiofrequency ablation and liver transplantation, its recurrence rate and prognosis remains unsatisfactory[2, 3]. Recently, immunotherapy has made great progress as a new treatment method in HCC. However, to date, this only benefited a subset of patients[4, 5]. The insufficient understanding of the tumor immunological microenvironment (TIME) may be the main reason for disappointing results. At the individual level, HCC has significant TIME heterogeneity, and the comprehensive understanding of the heterogeneity was crucial for clinical diagnosis, personalized treatment and prognosis prediction in HCC[6].

HCC is a typical inflammation-driven tumor, which is mainly derived from viral infections and liver fibrosis[7]. The transition from chronic hepatitis to HCC is accompanied by changes in local TIME [8]. Immune cells are the main components of TIME, and their number and status play a critical role in the progression of tumor development, invasion and metastasis. Previous researches have mainly focused on one or several immune cell types[9–12], which may bias the understanding of TIME due to the intensive cellular interaction between different cells. Hence, it is essential to be considered as a whole.

The rapid development of genomics and transcriptomics has made it possible to systematically explore the TIME heterogeneity in HCC. In the present study, we collected a total of 2,175 eligible samples from 15 cohorts, and combined with multi-omics data, hoping to explore different heterogeneous TIME phenotypes, further investigate the potential immune escape mechanisms of each TIME phenotype and what genomic alterations may lead to the formation of these different phenotypes. As a result, we successfully identified and validated three heterogeneous phenotypes based on the broad-spectrum immune cells in TIME. These three phenotypes exhibited different clinical outcomes, immune escape mechanisms and specific genomic alterations. In addition, the TIME index (TI) was developed to quantify TIME infiltration pattern, and it was a superior prognostic and immunotherapy predictor.

Data collecting and processing

The present workflow was shown in Fig. S1. For the discovery cohorts, the HCC microarray datasets were recruited from the Gene Expression Ominibus (GEO) database with the following criteria: (1) only from Affymetrix platform; (2) primary liver cancer; (3) untreated patients; (4) the number of patients was ≥50; (5) with more than 12,000 protein coding genes. Finally, 14 eligible datasets containing 1,821 patients were retrieved (Table S1). For the TCGA validation cohort, the TCGA-LIHC RNA-seq data was obtained from the UCSC Xena Portal. Please refer to Supplementary Materials and Methods for the data processing details. The corresponding clinical and sample information were obtained from the GEO and UCSC databases. For the TCGA-LIHC project, the somatic mutation data, copy number variation data, and DNA methylation data were obtained from the TCGA portal. In addition, we also downloaded the RNA-seq data and clinical information of the 32 other cancer types from the UCSC databases.

Integrated assessment of the TIME immune cell composition

In order to quantify the relative abundance of each immune cell population in TIME, we applied the single sample gene set enrichment analysis (ssGSEA) algorithm [13]. The gene sets for marking 24 immune cell types was recruited from Bindea et al. study [14]. In order to ensure the rationality and robustness of the ssGSEA results, two different algorithms were utilized to further validate: CIBERSORT [15] and MCP-counter [16]. The details were described in Supplementary Materials and Methods.

Identification and validation of the TIME phenotypes

We used the ConsensusClusterPlus package to determine the optimum number of clusters in the GEO cohort[17]. The results were further detected using the cumulative distribution function (CDF) curve, proportion of ambiguous clustering (PAC) score, and Nbclust [18]. To evaluate the reproducibility of the clusters generated from consensus clustering in the GEO cohort, the in-group proportion (IGP) statistical analysis was employed to further validate the existence of these clusters in the TCGA validation cohort [19]. The details were described in Supplementary Materials and Methods.

Gene set variation analysis (GSVA)

To further explore the potential biological function and progress variations of each phenotype, we conducted the GSVA analysis via GSVA package[20]. The gene sets, including the Hallmark and KEGG gene sets, were derived from the Molecular Signatures Database (MSigDB). The limma package was implemented to identify the significantly altered pathways between each phenotype and the others with the following threshold: log FC >0.2 and adjusted P-value <0.05. The resulting P-values from Benjamini-Hochberg (BH) multi-test correction were adjusted for multiple comparisons using the false discovery rate (FDR).

Analysis of immunogenomic features

We calculated or collected tumor mutation burden (TMB), SNV or Indel neoantigen load, aneuploidy scores (AS), homologous recombination defects (HRD) score, microsatellite instability (MSI), TCR or BCR diversity, cancer/testis-antigens (CTAs) level, antigen processing and presenting machinery scores (APS), and MHC-related molecules, in order to investigate the tumor immunogenicity of HCC. The details were described in Supplementary Materials and Methods.

Multi-omics profiling of immunomodulators

A total of 62 immunomodulators (including 12 MHC class I genes, 11 MHC class II genes, 27 checkpoint stimulator genes, and 12 checkpoint inhibitor genes) were recruited [21]. We investigated the multi-omics regulation landscape (including genes expression, somatic mutation, copy number variation (CNV), DNA methylation and miRNA expression) of 62 immunomodulators in three phenotypes. The Kruskal-Wallis test was performed for the mRNA expression, and Fisher’s test was performed for the somatic mutation and CNV of immunomodulators. The adjusted P-values were acquired using the BH multi-test correction. To survey the correlation between DNA methylation and gene expression of immunomodulators, each methylation site was matched to the corresponding gene. Most of the single genes had multiple methylation sites. In each phenotype, we assessed the Spearman’s correlation between each immunomodulator expression and all the corresponding methylation sites. Subsequently, we obtained a single correlation value for every gene by averaging the corresponding correlation coefficients. Next, we estimated the pattern that the miRNA modulated the immunomodulator expression. The significant inversely correlative pairs of miRNA and immunomodulator were included (Spearman correlation ≤ -0.2 and BH-corrected P <0.05) within each phenotype. Then, according to the predicted blinding targets for miRNA, these were curated from the miRDB database.

Driver mutation genes and mutation signatures

We utilized MutSigCV (version 1.41) to identify the significantly mutated genes (SMGs) for three phenotypes of TCGA-LIHC cohort[22], and a q-value of <0.05 was considered as the threshold. Subsequently, the MutationalPatterns R package was applied to extract the mutation signatures of each phenotype[23]. The mutational signatures can be extracted from mutation count matrix using non-negative matrix factorization (NMF). The optimal factorization rank, which was the number of mutational signatures, can be determined using the NMF package. After calculating the pairwise cosine similarity between the extracted mutation signatures and the 30 COSMIC signatures previously reported (http://cancer.sanger.ac.uk/cosmic/signatures), these extracted mutation signatures were then named based on the COSMIC signature.

Copy number variations

The TCGAbiolinks R package was used to download the CNV data based on the segment mean value (log2(copy-number/2)) obtained from TCGA database. The ABSOLUTE algorithm was implemented to estimate tumor ploidy for each sample [24]. To quantify the overall fraction of genomic alteration in three phenotypes, we calculated the fraction of genome alteration (FGA), fraction of genome gained (FGG), and fraction of genome lost (FGL). The FGA for a sample was defined as the ratio of the number of bases with CNVs to the number of all bases. The FGG or FGL considered only CNVs that were gained or lost. The GISTIC 2.0 was applied to define the recurrently amplified and deleted regions of each phenotype [25].

Methylation profiling

We downloaded the HumanMethylation450 array for HCC in TCGA. The global methylation level (GML) was estimated through averaged beta values of the specific probes, as described by Jung et al.[26]. The CD8+ T cell infiltrate status and proliferation score were both derived from Thorsson et al. study [21]. For each phenotype, we identified the epigenetically silenced genes (ESGs) using the following criteria: (1) excluding the CpG sites methylated in normal tissues (mean β-value of >0.2) or less than 10% of the tumor samples; (2) the DNA methylation data was divided into the methylation group and unmethylation group, according to the cutoff (β-value = 0.3); (3) for each probe, if the difference between the corresponding gene mean expression in the unmethylated group and that in the methylated group was >1.64 standard deviations of the unmethylated group, the probe would be labeled as epigenetically silenced; (4) when multiple probes were assigned to the same gene, the gene with more than half of the corresponding probes were labeled as epigenetically silenced, and identified as ESG.

TIME index

We applied the limma package to identify the differentially expressed genes (DEGs) between each phenotype and the others using the following thresholds: |log FC| >1.5 and adjusted P-value <0.05. Based on these DEGs, ssGSEA was performed to obtain the TIME index (TI) for each patient. Then, we assessed the prognostic value of TI in both the HCC and pancancer cohort. The performance of TI in predicting the response to immunotherapy was further evaluated. Six pre-treatment melanoma cohorts with the gene expression profile and immunotherapy information were collected (Supplementary Materials and Methods). The immunotherapy response prediction accuracy of the TI was compared with 11 other known biomarkers (CD274, PDCD1, CTLA4, CD8, TMB, T cell clonality, B cell clonality, TIDE, MSI Score, cytolytic activity (CYT)and APS; the details were described in Supplementary Materials and Methods). The area under the ROC curve was used as the quality metric of prediction. In addition, based on the nearest centroid method and Pearson’s correlation, we developed a pipeline to classify single patient into one of these three subtypes and calculated the TI (https://github.com/Zaoqu-Liu/TIME).

Immune cell infiltration patterns of TIME

We applied the ssGSEA method to assess the infiltration abundance of 24 immune cell types for 1,821 HCC samples. The correlation between these immune cells was presented in Fig. S1A. It was observed that several pairs strongly correlated with immune cells, such as T cell-cytotoxic cells, B cell-T cells and macrophage-immature DC cells. Subsequently, we performed a consensus cluster analysis, in which all HCC samples were initially grouped into different k (k = 2-9) clusters. The CDF curves of the consensus score and PAC value suggested that the optimal division was achieved when k = 3 (Fig. 2A-C). The same result was obtained from NbClust (Fig. S1B). The three clusters of samples were separated from each other on the two-dimensional principle component plot (Fig. 1D). Thus, based on the infiltration profiles of 24 immune cells in TIME, 1,821 HCC samples were finally classified into three TIME phenotypes (TIME-1 = 721, TIME-2 = 530, TIME-3 = 570). As shown in Fig. 1E-F, TIME-1 presented as an immune deficiency phenotype due to the lowest infiltration in almost all immune cells. On the contrary, it was found that TIME-3 had a significantly higher infiltration level in the majority of immune cells, especially adaptive immune cells (e.g. CD8+T cells and B cells), suggesting that TIME-3 was associated with immune activation and superior cytotoxic potential. TIME-2 was in an intermediate status of immune infiltration between TIME-1 and TIME-3, and was characterized by higher infiltration immunosuppressive cells that contain Treg and TH17. In addition, it was also observed that there was a richer infiltration in some innate immune cells that contain DC and NK cells in TIME-2. Two other different algorithms were further applied: CIBERSORT and MCP-counter. The results shared a consistent immune infiltration pattern with the ssGSEA method in HCC (Figs. S1C–S1D).

In order to ensure the reproducibility and robustness of the TIME phenotypes derived from the GEO cohort, we further conducted the IGP statistical method to validate the TIME phenotypes in the TCGA cohort. The three phenotypes were highly consistent between the discovery and validation cohorts, with the corresponding IGP values at 95.6%, 93.3% and 94.7%, respectively, and the three phenotypes were deemed to be of high-quality due to the statistically significance (all, P<0.001). The immune cells infiltration patterns in the TCGA cohort exhibited a very similar pattern of immune infiltration to the GEO cohort (Figs. S1E–S1F). Furthermore, the NbClust also indicated that the three clusters configuration was “optimal” in the TCGA cohort (Fig. S1G).

Specific functional pathways of each TIME phenotype

We further explored the specific functional status and biological mechanisms of each phenotype in the GEO cohort (Fig. 2G and Table S2). TIME-1 was prominently enriched in pathways, such as Myc targets, G2M checkpoint, and DNA repair. These pathways were remarkably associated with MKI67, and thereby with proliferation (Fig. 2H). Furthermore, it was observed that the immune-relevant pathways were significantly downregulated in TIME-1 (Table S2). Combined with the lack of immune cell infiltration, we inferred that TIME-1 may present an immune deficiency phenotype. On the contrary, TIME-3 enriched intensive pathways related to immune activation, and these pathways had a remarkably positive association with the immune score assessed through the ESTIMATE algorithm [27] (Fig. 2I), suggesting that TIME-3 may exhibit a state of immune activation. Notably, TIME-2 was significantly upregulated in metabolic-relevant pathways. It was observed that there was significantly negative correlation between the immune score and these specifically activated pathways in TIME-2 (Fig. 2I). Moreover, TIME-2 was rich in immunosuppressive cells (e.g. TH17 cell and Treg), which is known from the previous descriptions. Hence, it was concluded that TIME-2 may present as an immune-suppressed phenotype. The KEGG results was in accordance with the above (Fig. 2J and Table S3), and similar results were achieved in the TCGA cohort (Fig. S2A–S2B and Table S4-S5). Overall, we identified three TIME phenotypes in HCC showing significantly different immune cell infiltration and biological functions, respectively. TIME-1 was categorized as an immune-deficiency phenotype, characterized by immune cell depletion and proliferation; TIME-2 was categorized as an immune-suppressed phenotype, characterized by being in immunosuppressive states; TIME-3 was categorized as an immune-activated phenotype, characterized by abundant leukocytes infiltration and immune activation.

The clinical value of TIME phenotypes

We explored the prognostic value in the two independent cohorts (TCGA-LIHC and NCI cohort), which contained the complete overall survival (OS) and relapse free survival (RFS) information. The Kaplan-Meier analysis of both OS and RFS exhibited that HCC patients had an increasingly favorable prognosis from TIME-1 to TIME-3 (Figs. 3A-D). Furthermore, we assessed the sensitivity to sorafenib in TIME phenotypes by the pRRophetic package. TIME-1 was found to be more sensitive to sorafenib than the other phenotypes (Fig. 3E and Fig. S3A). Notably, TIME-1 exhibited the highest expression in sorafenib related target genes (Fig. 3F). Therefore, this suggested that patients in TIME-1 may benefit from sorafenib the most. As formerly mentioned, TIME-3 had higher levels of immune cell infiltration abundance (e.g. CD8+ T cells). Hence, we speculated that patients in TIME-3 might be more responsive to immunotherapy. The TIDE algorithm was applied to infer the response to immunotherapy. As respected, TIME-3 had a higher response rate than the other phenotypes (Fisher’s exact test: P=0.009) in the GEO cohort (Fig. 3G), and consistent results was found in the TCGA cohort (Fisher’s exact test: P=0.048) (Fig. S3B). We also utilized the submap algorithm to compare the similarity of the expression profiles between the three TIME phenotypes and 47 previous melanoma patients with detailed immunotherapeutic information, and revealed that patients in TIME-3 were more responsive to anti-PD1 treatment (Bonferroni corrected P=0.008)[28] (Fig. 3H). The submap analysis on the TCGA cohort also achieved similar results (Fig. S3C).

Potential extrinsic immune escape mechanism

To further research the regulatory mechanisms of the TIME phenotypes, we focused on the TCGA cohort, which possessed multiple omics data and comprehensive clinical data.

We firstly investigated the extrinsic immune escape mechanisms. Previous studies indicated that extrinsic immune escape may include three major aspects: absence of leukocytes, presence of immunosuppressive cells, and release of abundant immunosuppressive cytokines [29, 30]. As described above, TIME-1 was characterized by deficient immune cell infiltration and then lack of immune mediated elimination. TIME-2 was characterized by higher levels of immunosuppressive cells (e.g. TH17 cell and Treg; Figs. S4A–S4B), which indicated a role of immunosuppressive cells in immune escape. In addition, TIME-2 lacked immune active cells (e.g. CD8+T cells). Therefore, it was speculated that TIME-1 and TIME-2 probably reflect the absence of recruitment or activation of innate immune cells, inducing failure of adaptive anti-tumor immune responses. The low expression of molecules in TIME-1 and TIME-2, such as AIM2, TLR7 and TLR8, was potentially involved in priming of innate immunity, which further confirmed our speculation (Fig. 4A). TIME-3 was characterized by the presence of abundant innate and adaptive immune cells. In addition, TIME-3 had a higher expression of both immunostimulatory and immunoinhibitory cytokines, while these cytokines were all relatively low in TIME-1 and TIME-2 (Fig. S4C). These results implied that high concentrations of immunoinhibitory cytokines might contribute to the immune escape in TIME-3. It was noteworthy that the differential expression of cytokines in these three phenotypes could not be explained by the CNV and mutation frequency (all, P>0.05; Table S6). Overall, our analysis revealed that the extrinsic immune escape mechanisms of three phenotypes were lack of tumor-infiltrating leukocytes, increased immunosuppressive cells, and rich in immunoinhibitory cytokines, respectively.

Potential intrinsic immune escape mechanism

We further explored the potential intrinsic immune escape mechanisms in two major facets: tumor immunogenicity and the expression level of immune checkpoint molecules [31]. First, a series of elements associated with tumor immunogenicity were estimated: genomic instability degree, neoantigen burden, genomic viral content, CTA level, and tumor antigen presentation competence [32]. The first four elements were the main source of tumor-specific antigens. It was found that the genomic instability degree presented a decreasing trend from TIME-1 to TIME-3 (HRD, AS and MSI; all, P<0.05; Figs. 4B-D). Similarly, TIME-3 had a lower TMB, when compared to TIME-2 (P<0.05, Fig. 4E). In terms of genomic viral content, TIME-3 had more HBV read counts than TIME-1 (mean viral read counts, 44.365 vs. 20.907, P<0.05; Fig. S4D), in contrast to the HCV read counts (mean viral read counts, 12.931 vs. 16.111, P<0.1; Fig. S4E). Of note, the neoantigen burden was relatively lower in TIME-3, although the statistical difference among the three phenotypes in SNV or indel neoantigens was not reached (Figs. S4F–S4G). There was also not distinct variation among the CTAs overall expression of the three phenotypes (Fig. S4H). Overall, these above indicators had little difference in the TIME phenotypes. We further investigated the tumor antigen presentation capacity of the three phenotypes, and observed that TIME-3 had the highest APS and MHC-related molecules expression level, as opposed to TIME-1 (all, P<0.01; Fig. S4I and Fig. 4F), which was consistent with the CYT value and BCR/TCR diversity (all, P<0.05; Figs. 4G-I and Figs. S4J-S4K). This indicated that defective tumor antigen presentation capacity may be an intrinsic immune escape mechanism for TIME-1.

Subsequently, the genomic alterations of 62 immunomodulators were further summarized within the three TIME phenotypes (Fig. 4F). It was found that TIME-3 had higher costimulatory and coinhibitory molecules than the other phenotypes. This suggested that TIME-3 may overexpress the immune checkpoint molecules (such as CTLA4, CD274 and PDCD1; all BH-adjusted P<0.001) to evade the immune elimination after immune activation. All somatic mutations and CNVs did not significantly differ among the three phenotypes, and most of immunomodulators exhibited rare somatic mutations and CNVs, which indicated that the mutations and CNVs in the immunomodulators had little effect on TIME. Of note, DNA methylation negatively regulated many immunomodulators, such as CD27, CD226 and TNFSF8, implying epigenetic silencing. The associations were shown between miRNAs and immunomodulators for each TIME phenotype (BH corrected P-value <0.05; Fig. 4J), such as miR-17 negative correlated with CD274 and PDCD1LG2. It was also observed that the three phenotypes shared a common TNFSF4 negative regulator: miR-204. Compared with the mutation and CNV, methylation modification and miRNA sponges played leading roles in regulating the immunomodulators, indicating a new perspective for the development of immune checkpoint inhibitors.

Genomic alterations of the three TIME phenotypes

The investigators separately determined the SMGs among the three phenotypes using MutSigCV (Fig. 5A). All SMGs had mutation rates greater than 5% in three phenotypes. Among these three TIME phenotypes, the common SMGs (including TP53, CTNNB1 and ALB) had the top three significant MutSigCV q-value and frequent mutation rates, indicating that the mutation of TP53, CTNNB1 and ALB was broad in HCC. Additionally, the three phenotypes also displayed distinct SMGs, such as RB1, ACVR2A and CREB3L3 were SMGs of the three phenotypes, respectively. Besides, two newly identified SMGs, namely, BRD7 and RASA1, were classified as tumor suppressor genes, and these were associated with chromosome remodeling and cell proliferation[33, 34]. Based on the NMF, we isolated the mutation signatures of each phenotype. The age-related mutational processes (spontaneous deamination of 5-methylcytosine for signature 1 and unknow aetiology for signature 5) were prevalent in three phenotypes (Fig. 5B). TIME-1 had the least proportion than the other phenotypes (Fig. 5C). In addition, the three phenotypes also shared a common mutation profile (signature 22) associated with exposures to aristolochic acid. This may be associated with high-proportioned Asian patients in three phenotypes (Fig. S5A), and the aristolochic acid was mainly derived from herbal drugs of traditional Asian medicine[35]. Notably, signature 24, which represented the mutational pattern related with aflatoxin, was identified only in TIME-1, and possessed the maximum proportion (40.4%) (Fig. 5C). This implied that HCC patients in TIME-1 were more likely to be exposed to aflatoxin.

Tumor ploidy was estimated by ABSOLUTE, suggesting that a larger scale of HCC presented genome doubling, and that the doubling pattern was more frequent in TIME-1 compared with the other phenotypes (P=0.018, Fig. S5B). As shown in Fig. 4D-F, the FGA, FGG and FGL in TIME-1 were significantly higher than the other phenotypes, which might promote the cell proliferation and immune escape[36]. We further applied GISTIC 2.0 to delineate the significant focal copy number alterations of each phenotype (Fig. 5G and Table S7). CNVs that were recurrent in TIME-1 contained the focal amplification of 8q24.21 (MYC, ANXA13) and 13q34 (CDC16, TFDP1), and the focal loss of 14q22.1 (SAV1). Recurring focal arms CNVs in TIME-2 included the only amplification of 6p21.1 (VEGFA), and the loss of 13q13.3 (CCNV1). The genes on the focal loss arms of TIME-1 and TIME-2 were mainly associated with chemokines and cytokine through the GO annotation. Hence, the loss of these genes may contribute to the low immune infiltration of TIME-1 and TIME-2. TIME-3 exhibited the focal amplification that involved 8q24.12 (MTBP) and the focal deletion that involved 5q13.2 (TERT) and 10q23.31 (PTEN). These phenotype-special CNVs may play a crucial role in the biological features of the three phenotypes.

Furthermore, the combination of mutation and CNV data revealed the frequent alterations in different pathways (Fig. 5H). It was found that some genome alterations distributed evenly in these three phenotypes, such as TP53, CTNNB1 and ALB. In addition to these common alterations, we also observed the diverse alteration patterns in the pathways among three phenotypes. In TIME-1, cell cycle regulatory factor CDKN2A mutated in 12% of the cases. In TIME-2, VEGFA and its downstream genes PIC3KA and PTEN were frequently altered, and all of which were known to activate the PI3K pathway. APOB consumes an abundant cellular energy to facilitate very low-density lipoprotein (VLDL) secretion, and we found that APOB was frequently altered in 19% of the cases for TIME-2. Overall, the diverse genomic alteration preferences in three phenotypes might contribute to shape the TIME, and lead to differences in immune cell infiltration.

Methylation modification and regulation of the TIME phenotypes

As tumor cells divide, the loss of global methylation levels (GMLs) can result in chromosomal instability, and affect immune cell infiltration [37-39]. The GMLs in TIME-3 was significantly higher than the other phenotypes (P<0.001, Fig. 6A). Furthermore, there was a positive correlation between tumor-infiltrating CD8+ T cells and GMLs (Fig. 6B). In contrast, the decrease in GML might promote tumor cell proliferation (Fig. 6C). Subsequently, we identified 31, 25 and 39 ESGs in TIME-1, TIME-2 and TIME-3, respectively (Table S8 and Figs. 6D-E). It was noted that CPS1, FURIN and PHYHD1 were shared by the three phenotypes. As a liver-specific enzyme, CPS1 can facilitate cell division by instituting the pyrimidine synthesis pathway[40]. In TIME-1, FOXD4 was marked as epigenetically silenced in 85% of the cases, and its methylation silencing may lead to immune system dysfunction and tumor proliferation[41] (Fig. 6F). Tumor suppressors L1TD1 and PARP6, also the specific ESGs of TIME-1, and their methylation may promote tumor proliferation and low immune infiltration [42, 43] (Figs. 6G-H). SEMA3B belonged to the ESG of TIME-2, and its methylation can accelerate the progression of HCC[44]. FER1L is an another ESG of TIME-2, which can activate the PI3K/AKT pathway, further leading to the formation of the immunosuppressive microenvironment[45]. These results suggested that ESGs exhibit diversity among the three phenotypes, which may be involved in shaping the TIME, and these specific ESGs may also be potential therapeutic targets.

A robust prognostic and immunotherapy signature: TIME index

We identified 98 phenotype-related DEGs (Table S9). These genes, which significantly varied within the three phenotypes, possibly contributed to form the heterogenous TIME of HCC. In addition, many genes have been reported to be critical in immune response, such as CD27, CD8A, GZMA and IL7R [46-49]. The GO and KEGG annotation also displayed intensive immune phenotypes (Figs. S6A-S6B, Table S10-S11). Based on these DEGs, the ssGSEA was performed to obtain the TI of each patient. The TI presented a gradual increase from TIME1 to TIME3 (Figs. 7A-B). HCC patients with low TI were mainly distributed in TIME-1 (Fig. 7C). According to the optimal cut-off determined by the survminer package, we classified HCC patients into high and low TI groups. As expected, patients in the high TI group exhibited a tendency to better outcomes in the two independent cohort (Log-rank P<0.001, Figs. 7D-E). Besides, the multivariable Cox regression revealed that the TI was an independently prognostic factor in HCC (P=0.018, Table S12). We further assessed the TI of 10,121 patients that involved 33 differing types of cancers. Obvious TI diversity was observed in different cancers (Fig. 7F), which indicated that heterogeneous immune infiltration existed not only within the tumor, but also between tumors. The survival analysis for pancancer indicated that patients with a higher TI had a better prognosis (Log-rank P<0.001, Fig. 7G) and that the TI could independently affect the prognosis in the multivariable Cox regression (P=0.006, Table S12). Furthermore, it was found that 24 of 33 cancers presented a statistical significance in the Kaplan-Meier analysis (Fig. S7), and univariate Cox’s regression indicated that the TI was a protective factor in many different tumors (Fig. 7H).

We further explored the predictive ability of the TI for immunotherapy. In order to determine the predictive power, we also computed the response prediction of 11 other known biomarkers. The area under the ROC curve (AUC) was used as the quality metric of prediction. We found that the TI exhibited robust predictions in the six cohorts (Fig. 7I). Especially in the Riaz et al. cohort, the TI reached the highest prediction accuracy with AUC = 0.955. It was also observed that although the AUC value of TI in the Liu et al. cohort was less than 0.7, this was still the highest, when compared to the other markers. In addition, the TI displayed the highest mean AUC value compared with the other biomarkers (Fig. 7J). Hence, the present work strongly suggests that the TI was a potential and robust biomarker for the prognosis and clinical response assessment of immunotherapy.

In addition, to advance clinical application, we developed a R package termed “TIME” in the GitHub website (https://github.com/Zaoqu-Liu/TIME). The pipeline could classify single patient into one of these three subtypes and calculated the TI.

The hepatocellular carcinoma (HCC) ecosystem is diverse, complex and dynamic in nature, and is mainly composed of tumor cells and immune cells[50]. Immune cells are the main components of TIME, and their number and status play a crucial role in the progression of tumor development, invasion and metastasis. To the best of our knowledge, the present study is the first to systematically investigate the heterogeneity of TIME from the dimension of broad-spectrum immune cells, and comprehensively explore the potential immune escape mechanisms and specific genomic alterations of different TIME phenotypes. In addition, the TI was proposed to quantify TIME infiltration pattern, and it was also a superior prognostic and immunotherapy predictor. These results can enhance the understanding of TIME, and guide more effective personalized immunotherapies.

As described, TIME-1 was characterized by immune cell depletion and proliferation, corresponding to the immune-deficiency phenotype, TIME-2 was characterized by being in immunosuppressive states, corresponding to the immune-suppressed phenotype, and TIME-3 was characterized by abundant immune cell infiltration and immune activation, corresponding to immune-activated phenotype. We underlined the potential immune escape mechanism of each phenotype: lack of leukocytes and defective tumor antigen presentation capacity in TIME-1, increased immunosuppressive cells in TIME-2, and rich in immunoinhibitory molecules in TIME-3. The distinct immune escape mechanism among TIME phenotypes could provide strategies for improving the efficacy of HCC immunotherapy. In addition, it was also found that methylation modification and miRNA sponges may play leading roles in regulating the immunomodulator expression, when compared to mutations and CNVs, suggesting that the development of ICIs should consider the methylation modification and miRNA regulation of immunomodulators.

The TIME phenotypes also exhibited a significant heterogeneity at the genomic level, which may drive the formation of different phenotypes. TIME-1 was dominated by the alterations of TP53, CDKN2A, CTNNB1 and AXIN1. CDKN2A could inhibit the tumor-promoting behavior of CDK4/6, and given that the CDKN2A alteration is frequent in HCC, CDK4/6 inhibitors are presently being tested in advanced HCC[51]. Thus, it was suspected that TIME-1 may be more sensitive to CDK4/6 inhibitors. Previous studies showed that the mutations of CTNNB1 and AXIN1 might be characteristic of immune exclusion, and represent the biomarkers of innate resistance to immunotherapy[52, 53]. Besides, the FOXD4 methylation was also associated with the immune dysfunction and cell proliferation. These genomic events might contribute to the immunodeficiency and proliferation of TIME-1. TIME-2 was characterized by significant alteration patterns in the PI3K pathway, such as PIK3CA, VEGFA and PTEN. The mutation of PIK3CA may serve as a reliable biomarker for relatively poor response to immunotherapies, such as PD-L1 antibodies [54]. VEGFA, as a tumor angiogenic factor, also plays a pivotal role in the formation of the tumor immunosuppressive microenvironment. Its alterations induce the proliferation of TH17 and Treg, and inhibit CD8 + T cell function, resulting in immune escape[55–57], which in line with the immune cell infiltration characteristic of TIME-2. Of note, TIME-2 might be more sensitive to the VEGFA-targeting monoclonal antibody Bevacizumab due to its significant VEGFA amplification[58]. ARID1A was the SMG of TIME-3, and its inactivating mutations can lead to remarkable increases in CD8 and PD-L1, and tumors with ARID1A deficiency were also more sensitive to PD-L1 antibodies[59]. Overall, these results suggested that distinct genomic alterations might not only lead to different immune cell infiltration and functional status, but also explain the potential reasons for the sensitivity or resistance of different phenotypes to immunotherapy, which provide references for the precise treatment of HCC.

The TIME phenotypes have significant clinical value. Consistent with the immune infiltration of three phenotypes, TIME-1 indicated the worst prognosis, while TIME-3 had the most favorable prognosis. TIME-1 was predicted to be most sensitive to sorafenib, which was consistent with higher expression of drug targets. Unsurprisingly, in line with the higher level of immune cells infiltration and immune checkpoint molecule expression, TIME-3 exhibited a superior response to immunotherapy. In addition, we proposed a scoring scheme to quantify the TIME infiltration pattern termed TIME index (TI). The TI was not only an independent prognostic biomarker for both HCC and pancancer, but also performed well in predicting the response to immunotherapy. Hence, the TI could guide clinical management and personalized immunotherapy of HCC. In addition, for the purpose of facilitating clinical application, we developed a pipeline to classify single patient into one of these three subtypes and calculated the TI.

The present study also had some limitations. First, we only considered the inter-individual heterogeneity due to the lack of data, but did not consider the intra-tumor heterogeneity, which is common in multifocal HCC[60]. Second, although machine learning algorithms were applied to predict the sensitivity of the TIME phenotypes to sorafenib and immunotherapy, further clinical validation is need. Finally, potential genomic drivers require further functional verification. Clinical studies and related experiments are ongoing in our hospitals and laboratories.

In summary, our research revealed the three heterogeneous TIME phenotypes with different clinical outcomes, immune escape mechanisms, and genomic alterations in HCC, which could present strategies for improving the efficacy of immunotherapy. The TI as a novel prognostic and immunotherapeutic signature that could guide clinical management and personalized immunotherapy.

TIME: Tumor immunological microenvironment; HCC: Hepatocellular carcinoma; TI: TIME index; GEO: Gene Expression Ominibus; ssGSEA: Single sample gene set enrichment analysis; CDF: Cumulative distribution function; PAC: Proportion of ambiguous clustering; IGP: In-group proportion; GSVA: Gene set variation analysis; TMB: Tumor mutation burden; AS: Aneuploidy scores; HRD: Homologous recombination defects; MSI: Microsatellite instability; CTAs: Cancer/testis-antigens; APS: Antigen processing and presenting machinery scores; CNV: Copy number variation; SMGs: Significantly mutated genes; NMF: Non-negative matrix factorization; FGA: Fraction of genome alteration; FGG: Fraction of genome gained; FGL: Fraction of genome lost; GML: Global methylation level; ESGs: Epigenetically silenced genes; DEGs: Differentially expressed genes; OS: Overall survival; RFS: Relapse free survival; CYT: Cytolytic activity; VLDL: Very low-density lipoprotein; AUC: Area under the ROC curve.

Ethics approval and consent to participate

The patient data in this work were acquired from the publicly available datasets whose informed consent of patients were complete.

Consent for publication

All Authors have seen and approved the manuscript and consent publication.

Availability of data and materials

All data used in this work can be acquired from the Gene-Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) and the GDC portal (https://portal.gdc.cancer.gov/).

Competing interests

The authors declare that they have no competing interests.

Funding

Not applicable.

Authors’ contributions

ZQL, XWH and ZQS designed this work. ZQL, DCJ, YYZ and CCS integrated and analyzed the data. ZQL, CCS and YYZ wrote this manuscript. XLZ, ZQL, KHX and ZQS edited and revised the manuscript. All authors approved this manuscript.

Acknowledgements

Not applicable.

Duran SR, Jaquiss RDB. Hepatocellular Carcinoma. N Engl J Med. 2019; 381: e2.
Iizuka N, Hamamoto Y, Oka M. Predicting individual outcomes in hepatocellular carcinoma. Lancet. 2004; 364: 1837-9.
Iizuka N, Oka M, Yamada-Okabe H, Nishida M, Maeda Y, Mori N, et al. Oligonucleotide microarray for prediction of early intrahepatic recurrence of hepatocellular carcinoma after curative resection. Lancet. 2003; 361: 923-9.
Finn RS, Qin S, Ikeda M, Galle PR, Ducreux M, Kim TY, et al. Atezolizumab plus Bevacizumab in Unresectable Hepatocellular Carcinoma. N Engl J Med. 2020; 382: 1894-905.
Hou J, Zhang H, Sun B, Karin M. The immunobiology of hepatocellular carcinoma in humans and mice: Basic concepts and therapeutic implications. J Hepatol. 2020; 72: 167-82.
Zhang Q, Lou Y, Yang J, Wang J, Feng J, Zhao Y, et al. Integrated multiomic analysis reveals comprehensive tumour heterogeneity and novel immunophenotypic classification in hepatocellular carcinomas. Gut. 2019; 68: 2019-31.
Ghouri YA, Mian I, Rowe JH. Review of hepatocellular carcinoma: Epidemiology, etiology, and carcinogenesis. J Carcinog. 2017; 16: 1.
Ringelhan M, Pfister D, O'Connor T, Pikarsky E, Heikenwalder M. The immunology of hepatocellular carcinoma. Nat Immunol. 2018; 19: 222-32.
Zheng C, Zheng L, Yoo JK, Guo H, Zhang Y, Guo X, et al. Landscape of Infiltrating T Cells in Liver Cancer Revealed by Single-Cell Sequencing. Cell. 2017; 169: 1342-56 e16.
Garnelo M, Tan A, Her Z, Yeong J, Lim CJ, Chen J, et al. Interaction between tumour-infiltrating B cells and T cells controls the progression of hepatocellular carcinoma. Gut. 2017; 66: 342-51.
Ouyang FZ, Wu RQ, Wei Y, Liu RX, Yang D, Xiao X, et al. Dendritic cell-elicited B-cell activation fosters immune privilege via IL-10 signals in hepatocellular carcinoma. Nat Commun. 2016; 7: 13453.
Peng H, Wisse E, Tian Z. Liver natural killer cells: subsets and roles in liver immunity. Cell Mol Immunol. 2016; 13: 328-36.
Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009; 462: 108-12.
Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013; 39: 782-95.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12: 453-7.
Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016; 17: 218-.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26: 1572-3.
Șenbabaoğlu Y, Michailidis G, Li JZ. Critical limitations of consensus clustering in class discovery. Sci Rep-Uk. 2014; 4: 6207-.
Kapp AV, Tibshirani R. Are clusters found in one dataset present in another dataset? Biostatistics. 2007; 8: 9-31.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14: 7-.
Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The Immune Landscape of Cancer. Immunity. 2018; 48: 812-30 e14.
Lawrence MS, Stojanov P, Polak P, Kryukov GV, Cibulskis K, Sivachenko A, et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013; 499: 214-8.
Blokzijl F, Janssen R, van Boxtel R, Cuppen E. MutationalPatterns: comprehensive genome-wide analysis of mutational processes. Genome Med. 2018; 10: 33.
SL C, K C, E H, A M, H S, T Z, et al. Absolute quantification of somatic DNA alterations in human cancer. Nature biotechnology. 2012; 30: 413-21.
CH M, SE S, B H, ML M, R B, G G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome biology. 2011; 12: R41.
Jung H, Kim HS, Kim JY, Sun J-M, Ahn JS, Ahn M-J, et al. DNA methylation loss promotes immune evasion of tumours with high mutation and copy number load. Nature Communications. 2019; 10: 4278.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4: 2612.
Roh W, Chen PL, Reuben A, Spencer CN, Prieto PA, Miller JP, et al. Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci Transl Med. 2017; 9.
Schreiber RD, Old LJ, Smyth MJ. Cancer immunoediting: integrating immunity's roles in cancer suppression and promotion. Science. 2011; 331: 1565-70.
Beatty GL, Gladney WL. Immune escape mechanisms as a guide for cancer immunotherapy. Clin Cancer Res. 2015; 21: 687-92.
Schreiber RD, Old LJ, Smyth MJ. Cancer immunoediting: integrating immunity’s roles in cancer suppression and promotion. Science. 2011; 331: 1565-70.
Schumacher TN, Schreiber RD. Neoantigens in cancer immunotherapy. Science. 2015; 348: 69-74.
Drost J, Mantovani F, Tocco F, Elkon R, Comel A, Holstege H, et al. BRD7 is a candidate tumour suppressor gene required for p53 function. Nat Cell Biol. 2010; 12: 380-9.
Maertens O, Cichowski K. An expanding role for RAS GTPase activating proteins (RAS GAPs) in cancer. Adv Biol Regul. 2014; 55: 1-14.
Forbes SA, Beare D, Boutselakis H, Bamford S, Bindal N, Tate J, et al. COSMIC: somatic cancer genetics at high-resolution. Nucleic Acids Res. 2017; 45: D777-D83.
T D, H U, EC W, SJ E. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science (New York, NY). 2017; 355.
Eden A, Gaudet F, Waghmare A, Jaenisch R. Chromosomal instability and tumors promoted by DNA hypomethylation. Science. 2003; 300: 455.
Davoli T, Uno H, Wooten EC, Elledge SJ. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science. 2017; 355.
Zhou W, Dinh HQ, Ramjan Z, Weisenberger DJ, Nicolet CM, Shen H, et al. DNA methylation loss in late-replicating domains is linked to mitotic cell division. Nat Genet. 2018; 50: 591-602.
Liu H, Dong H, Robertson K, Liu C. DNA Methylation Suppresses Expression of the Urea Cycle Enzyme Carbamoyl Phosphate Synthetase 1 (CPS1) in Human Hepatocellular Carcinoma. 2011; 178: 652-61.
Jackson BC, Carpenter C, Nebert DW, Vasiliou V. Update of human and mouse forkhead box (FOX) gene families. Hum Genomics. 2010; 4: 345-52.
Altenberger C, Heller G, Ziegler B, Tomasich E, Marhold M, Topakian T, et al. SPAG6 and L1TD1 are transcriptionally regulated by DNA methylation in non-small cell lung cancers. Mol Cancer. 2017; 16: 1.
Honda S, Minato M, Suzuki H, Fujiyoshi M, Miyagi H, Haruta M, et al. Clinical prognostic value of DNA methylation in hepatoblastoma: Four novel tumor suppressor candidates. Cancer Sci. 2016; 107: 812-9.
Tischoff I, Markwarth A, Witzigmann H, Uhlmann D, Hauss J, Mirmohammadsadegh A, et al. Allele loss and epigenetic inactivation of 3p21.3 in malignant liver tumors. Int J Cancer. 2005; 115: 684-9.
Wang X, Dong K, Jin Q, Ma Y, Yin S, Wang S. Upregulation of lncRNA FER1L4 suppresses the proliferation and migration of the hepatocellular carcinoma via regulating PI3K/AKT signal pathway. J Cell Biochem. 2019; 120: 6781-8.
Denoeud J, Moser M. Role of CD27/CD70 pathway of activation in immunity and tolerance. Journal of leukocyte biology. 2011; 89: 195-203.
Suetake H, Araki K, Akatsu K, Somamoto T, Dijkstra JM, Yoshiura Y, et al. Genomic organization and expression of CD8alpha and CD8beta genes in fugu Takifugu rubripes. Fish & shellfish immunology. 2007; 23: 1107-18.
Shimizu K, Yamasaki S, Sakurai M, Yumoto N, Ikeda M, Mishima-Tsumagari C, et al. Granzyme A Stimulates pDCs to Promote Adaptive Immunity via Induction of Type I IFN. Frontiers in immunology. 2019; 10: 1450.
Guo Z, Johnson V, Barrera J, Porras M, Hinojosa D, Hernández I, et al. Targeting cytochrome P450-dependent cancer cell mitochondria: cancer associated CYPs and where to find them. Cancer metastasis reviews. 2018; 37: 409-23.
Liu J, Dang H, Wang XW. The significance of intertumor and intratumor heterogeneity in liver cancer. Exp Mol Med. 2018; 50: e416.
Knudsen ES, Gopal P, Singal AG. The changing landscape of hepatocellular carcinoma: etiology, genetics, and therapy. Am J Pathol. 2014; 184: 574-83.
Llovet JM, Montal R, Sia D, Finn RS. Molecular therapies and precision medicine for hepatocellular carcinoma. Nat Rev Clin Oncol. 2018; 15: 599-616.
Pinyol R, Sia D, Llovet JM. Immune Exclusion-Wnt/CTNNB1 Class Predicts Resistance to Immunotherapies in HCC. Clin Cancer Res. 2019; 25: 2021-3.
Choi M, Kadara H, Zhang J, Parra ER, Rodriguez-Canales J, Gaffney SG, et al. Mutation profiles in early-stage lung squamous cell carcinoma with clinical follow-up and correlation with markers of immune function. Ann Oncol. 2017; 28: 83-9.
Suzuki T, Hirakawa S, Shimauchi T, Ito T, Sakabe J, Detmar M, et al. VEGF-A promotes IL-17A-producing gammadelta T cell accumulation in mouse skin and serves as a chemotactic factor for plasmacytoid dendritic cells. J Dermatol Sci. 2014; 74: 116-24.
Mandal R, Senbabaoglu Y, Desrichard A, Havel JJ, Dalin MG, Riaz N, et al. The head and neck cancer immune landscape and its immunotherapeutic implications. JCI Insight. 2016; 1: e89829.
Voron T, Colussi O, Marcheteau E, Pernot S, Nizard M, Pointet AL, et al. VEGF-A modulates expression of inhibitory checkpoints on CD8+ T cells in tumors. J Exp Med. 2015; 212: 139-48.
Garcia J, Hurwitz HI, Sandler AB, Miles D, Coleman RL, Deurloo R, et al. Bevacizumab (Avastin(R)) in cancer treatment: A review of 15 years of clinical experience and future outlook. Cancer Treat Rev. 2020; 86: 102017.
Shen J, Ju Z, Zhao W, Wang L, Peng Y, Ge Z, et al. ARID1A deficiency promotes mutability and potentiates therapeutic antitumor immunity unleashed by immune checkpoint blockade. Nat Med. 2018; 24: 556-62.
Feo F, Pascale RM. Multifocal hepatocellular carcinoma: intrahepatic metastasis or multicentric carcinogenesis? Ann Transl Med. 2015; 3: 4.

SupplementTable.xlsx
SupplementaryMaterialsandMethods.docx
SupplementaryFig.S1.jpg
The immune infiltration pattern of the TIME phenotypes. A, Association between 24 immune cell subsets. B, Recommended number of clusters using 26 criteria of Nbclust package in the GEO cohort. C, The immune cells infiltration pattern assessed by CIBERSORT algorithm. D, The immune cells infiltration pattern assessed by MCP-counter algorithm. E, The infiltration abundance of 24 immune cell subsets evaluated by ssGSEA algorithm for three TIME phenotypes in TCGA cohort. Survival status, age, gender, vascular invasion, histology grade, AJCC stage, the TIME phenotypes are shown as patient annotations. F, The differences of 24 immune cell subsets infiltration among the three TIME phenotypes in the TCGA cohort. G, Recommended number of clusters using 26 criteria of Nbclust package in the TCGA cohort. For the boxplot, the asterisks represented the statistical p value (*P < 0.05, **P < 0.01, *** P < 0.001, **** P < 0.0001).
SupplementaryFig.S2.jpg
The specific functional status and biological characteristics of each TIME phenotype in the TCGA cohort. A, The activation states of Hallmark pathways of distinct TIME phenotypes in the TCGA cohort. B, The activation states of KEGG pathways of distinct TIME phenotypes in the TCGA cohort.
SupplementaryFig.S3.jpg
Assessment of chemotherapy and immunotherapy in the TCGA Cohort. A, Distribution of the estimated IC50 of sorafenib among three TIME phenotypes in the TCGA cohort. B, The distribution of the immunotherapy response results predicted by TIDE algorithm among three TIME phenotypes in the TCGA cohort. C, Submap analysis of the TCGA cohort and 47 previous melanoma patients with detailed immunotherapeutic information. For the boxplot, the asterisks represented the statistical p value (*P < 0.05, **P < 0.01, *** P < 0.001, **** P < 0.0001).
SupplementaryFig.S4.jpg
Potential immune escape mechanisms of each phenotype. A-B, The difference of TH17 (A) and Treg (B) among the three TIME phenotypes. C, The mRNA expression of chemokines, interleukins, interferons, and other important cytokines and their receptors for each TIME phenotype. D-K, The distribution of HBV read counts (D), HCV read counts (E), SNV neoantigens (F), indel neoantigens (G), CTA score (H), APS value (I), BCR-Richness diversity (J) and TCR-Richness diversity (K) in three TIME phenotypes. For the boxplot, the asterisks represented the statistical p value (nsP > 0.05, *P < 0.05, **P < 0.01, *** P < 0.001, **** P < 0.0001). For the heatmap, the asterisks represented the statistical p value (*P < 0.05, **P < 0.01, *** P < 0.001, **** P < 0.0001).
SupplementaryFig.S5.jpg
The distribution of race and ploidy. A, The distribution of race in the TCGA cohort. B, The distribution of tumor ploidy in three TIME phenotypes.
SupplementaryFigs.S6.jpg
GO and KEGG enrichment analysis of differentially expressed genes (DEGs).
SupplementaryFig.S7.jpg
Kaplan–Meier curves of OS for the TIME index (TI) in the 33 cancer types.

Download PDF

Journal Publication

published 06 Jan, 2021

Read the published version in Journal of Translational Medicine →

Editorial decision: Major revision
22 Dec, 2020
Review #2 received at journal
21 Dec, 2020
Reviewer #3 agreed at journal
09 Dec, 2020
Review #1 received at journal
23 Nov, 2020
Reviewer #2 agreed at journal
09 Nov, 2020
Reviewers invited by journal
09 Nov, 2020
Reviewer #1 agreed at journal
09 Nov, 2020
Editor assigned by journal
04 Nov, 2020
Editor invited by journal
04 Nov, 2020
Submission checks completed at journal
29 Oct, 2020
First submitted to journal
26 Oct, 2020

You are reading this latest preprint version

A Novel Immune Classification Reveals Distinct Immune Escape Mechanism and Genomic Alterations: Implications for Immunotherapy in Hepatocellular Carcinoma

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Methods

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1