Research Paper Volume 12, Issue 19 pp 19641—19659
Identification of immune landscape signatures associated with clinical and prognostic features of hepatocellular carcinoma
- 1 Department of Radiation Oncology, Nanfang Hospital, Southern Medical University, Guangzhou 510515, China
- 2 Division of Hepatobiliopancreatic Surgery, Department of General Surgery, Nanfang Hospital, Southern Medical University, Guangzhou 510515, China
- 3 The First School of Clinical Medicine, Southern Medical University, Guangzhou 510515, China
- 4 Medical Imaging Specialty, the First School of Clinical Medicine, Southern Medical University, Guangzhou 510515, China
- 5 Department of Pathology, Nanfang Hospital, Southern Medical University, Guangzhou 510515, China
- 6 Clinical Medicine Specialty, the First School of Clinical Medicine, Southern Medical University, Guangzhou, 510515, China
Received: February 21, 2020 Accepted: August 14, 2020 Published: October 13, 2020
https://doi.org/10.18632/aging.103977How to Cite
Copyright: © 2020 Yan et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Abstract
While cancer immunotherapy has been remarkably successful in some malignancies, some cancers derive limited benefit from current immunotherapies. Here, we combined immune landscape signatures with hepatocellular carcinoma clinical and prognostic features to classify them into distinct subtypes. The immunogenomic profiles, stromal cell features and immune cell composition of the subtypes were then systematically analyzed. Two independent prognostic indexes were established based on 6 immune-related genes and 17 differentially expressed genes associated with stromal cell content. These indexes were significantly correlated with tumor mutation burden, deficient DNA mismatch repair and microsatellite instability. In addition, tumor-infiltrating lymphocytes, including activated NK cells, resting memory CD4 T-cells, eosinophils, and activated mast cells were significantly correlated with hepatocellular carcinoma survival. In conclusion, we have comprehensively described the immune landscape signatures and identified prognostic immune-associated biomarkers of hepatocellular carcinoma. Our findings highlight potential novel avenues for improving responses to immunotherapy.
Introduction
Hepatocellular carcinoma (HCC) is ranked as the 3rd leading cause of cancer deaths globally [1]. Most HCCs are associated with liver cirrhosis arising from chronic infection of hepatitis B and C viruses. Liver cirrhosis is also caused by alcohol consumption and fatty liver disease [2]. Treatment strategies for HCC are dependent on disease stage. Surgery is the standard treatment for early stage HCC and has a 70% 5-year survival rate. Where surgery or liver transplantation is unfeasible, loco-regional therapies, including radiotherapy, radiofrequency, thermal and non-thermal ablation, and trans-arterial chemoembolization may be recommended as second line therapies. However, the 3-5-year survival rates of such approaches are highly variable [3]. For advanced unresectable HCC, the recommended treatments include regorafenib and sorafenib (tyrosine kinase inhibitors (TKIs)) and the vascular endothelial growth factor receptors 1 through 3 (VEGFR1-3) inhibitor, lenvatinib. However, these therapies offer very limited survival benefit [4].
Cancer immunotherapy has emerged as one of the most promising cancer treatments for various types of advanced solid tumors, including HCC [5]. Several immune-based therapies such as peptide vaccines against HCC antigens, natural killer cell therapies, and chimeric antigen receptor-engineered T cell therapies are under investigation in phase I/II trials [5]. Among them, immune checkpoint inhibitors have shown good prospects in the treatment of advanced HCC as evidenced by durable objective response rates (ORRs) and acceptable safety profiles in phase I/II trials [6].
While immunotherapy is remarkably successful in many HCC cases, it remains unclear why a subset of HCC patients fails to respond to it. Intratumoral genetic heterogeneity is a well-known feature of cancer [7]. Tumor mutation burden (TMB), neoantigen load, programmed cell death ligand 1 (PD-L1) levels, impaired DNA mismatch repair (MMR), and microsatellite instability (MSI) have been found to affect response of cancer cells to immunotherapy [8, 9]. Additionally, the highly heterogeneous immune microenvironment has been associated with various aspects of cancer [10]. For instance, non-tumor cell types within or around tumors, including stromal and immune cells, also influence cancer progression [11]. Thus, development of potent immunotherapy requires a comprehensive understanding of the heterogeneity of tumor and tumor immune microenvironment (TIME).
In the present study, HCC tissues were subclassified on the basis of immunogenomic profiles, stromal cell features, and immune cell composition. We systematically analyzed the molecular features of each sub-class, such as gene ontology, genes, chemotactic factors, regulatory pathways and networks, and integrated them with various clinicopathological features and patient clinical outcomes. Our findings highlight the potential clinical utility of individualized immune signatures in prognostic stratification and personalized immunotherapy for HCC patients.
Results
Immunogenomic profiling identifies 2 HCC subtypes
Enrichment level and activity of several immune cells, pathways or functions in HCC were analyzed using single sample gene-set enrichment analysis (ssGSEA) score based on 29 immune-associated gene sets [12, 13]. According to the ssGSEA score of these gene sets, we hierarchically clustered 369 HCC samples from TCGA datasets. Consequently, 2 distinct clusters, termed Immunity_L (Immunity Low) and Immunity_H (Immunity High) were identified (Figure 1A). Lower immune scores were seen in Immunity_L than in Immunity_H (Figure 1B). Notably, analysis of tumor purity and stromal score revealed opposite results, with tumor purity being low in Immunity_H and high in Immunity_L, while stromal score was high in Immunity_H and low in Immunity_L (Figure 1C, 1D). This suggests that Immunity_L group harbors more tumor cells whereas Immunity_H harbors more stromal and immune cells.
Figure 1. Characterization of 2 HCC subtypes based on immunogenomic profiling. (A) Heatmap of normalized ssGSEA scores for tumor purity, ESTIMATE score, Immune Score, and Stromal Score. Tumor samples were grouped into 2 immune classes, Immunity_H and Immunity_L, based on 29 immune-associated gene sets. (B–D) Violin plot analysis comparing the Immune Score (B), Stromal Score (C) and tumor purity (D) between Immunity_H and Immunity_L (Mann-Whitney U test). (E) Box plot comparison of relative HLA gene expression between HCC subtypes (T-test). (F) Box plots comparing expression levels of various immunologic activity related chemotactic factors between HCC subtypes. (G) Comparison of PD-L1 expression in Immunity_H and Immunity_L (T-test). (H) Comparison of 3-year survival between HCC subtypes in TCGA and GSE14520 datasets (log-rank test).
It was seen that most HLA (human leukocyte antigen) genes were upregulated in Immunity_H than in Immunity_L (Figure 1E). Levels of various immunologic chemotactic factors were also markedly elevated in Immunity_H, including T-cell, monocyte, macrophage, mast cell and eosinophil chemotactic factors (Figure 1F) [14]. Evaluation of PD-L1 expression in both HCC subtypes showed that PD-L1 levels were markedly higher in Immunity_H than in Immunity_L (Figure 1G). This implies that the Immunity_H subtype may show a stronger response to anti-PD-L1 immunotherapy given that PD-L1 expression positively associates with immunotherapeutic responsiveness [15].
Survival analyses of 2 public datasets showed that the 2 HCC subtypes were associated with different clinical outcomes. The Immunity_H subtype corresponded with significantly better 3-year survival relative to Immunity_L (Figure 1H), which is consistent with previous findings in which HCC associated with higher immune activity exhibited better clinical outcomes [16]. These data highlight the value of subtyping as the basis of immunogenomic profiling.
Changes in TIME and prognosis based on differential immune-related genes
Identification of HCC subtype-specific gene expression profiles
We then used EdgeR to identify HCC subtype-specific gene expression profiles [17]. A total of 729 genes were differentially expressed between Immunity_H and Immunity_L subtypes. Among these, 706 were upregulated and 23 were downregulated (Figure 2A, 2B). GO term analysis revealed that MHC-class II protein complex, natural killer cell chemotaxis and positive regulation of interrleukin-2 biosynthetic process were the most significantly enriched cellular components and biological process in Immunity_H (Figure 2C, 2D). KEGG pathway analysis identified allograft rejection, antigen presentation and processing to be most highly enriched in Immunity_H (Figure 2E, F). These genetic alterations affect the immune system modulations [18–21].
Figure 2. HCC subtype-specific gene expression profiles. (A, B) Heatmap and volcano plot of differentially expressed genes between Immunity_H and Immunity_L. Red, green and black dots indicate upregulated, downregulated and unchanged genes. It also applies to the following figures. (C–F) Significantly enriched GO terms (C, D) and KEGG pathways (E, F) of HCC subtypes.
Identification of differentially expressed immune-related genes
From the above set of genes, we identified 181 differentially expressed immune-related genes (IRGs) between Immunity_H and Immunity_L, of which 180 were upregulated and 1 was downregulated (Figure 3A, 3B). To identify the IRGs actively participating in HCC tumorigenesis and progression, we selected 22 differentially expressed IRGs that showed marked correlation with clinical outcomes. A forest plot of hazard ratios revealed 3 of these IRGs as protective factors, and the remaining 19 as predictors of poor prognosis (Figure 3C). To elucidate the molecular mechanisms underlying the clinical effects of the 22 IRGs, we evaluated the expression patterns of 318 transcription factors (TFs). Consequently, 13 TFs were found to be substantially higher in Immunity_H than in Immunity_L (Figure 3D, 3E). A regulatory network was established using the IRGs and TFs. A TF-based regulatory schematic was used to show the regulatory relationships between 7 of the 13 TFs and 20 of the 22 prognosis-associated IRGs (Figure 3F).
Figure 3. HCC subtype-specific differentially expressed IRGs. (A, B) Heatmap and volcano plot of the differentially expressed IRGs between Immunity_H and Immunity_L. (C) Forest plot of hazard ratios and corresponding 95% confidence intervals were estimated from univariate Cox’s regression analyses. Variables significantly associated with a good and poor OS are shown in green and red, respectively. (D, E) Volcano plot and heatmap of the differentially expressed TF genes between Immunity_H and Immunity_L. (F) Combinatorial TFs-IRGs regulation networks. Red, blue, and green, indicate high risk genes, low risk genes, and TF genes. It also applies to the following figures.
Evaluation of clinical outcomes based on the prognostic IRG panel
Multivariate Cox regression analysis was performed to assess the risk score for the prognosis-associated IRGs and constructed a prognostic signature for stratification of HCC patients into 2 groups based on the clinical outcomes. Disease-free survival, overall survival, cancer specific survival and progression-free survival were higher in the low risk group than in the high risk group (Figure 4A). The area under the curve (AUC) of the survival-dependent receiver operating characteristic (ROC) curve was 0.761, which was higher than that of clinicopathologic factors, indicating a high prognostic power for this IRG panel (Figure 4B). Heatmaps were established to visualize the expression profile of the 6 genes included in this panel for the low and high risk groups (Figure 4C). The risk score was calculated as follows: [Expression level of IL18RAP * (-2.5273)] + [Expression level of FCER1G * 0.0038] + [Expression level of CXCL11 * 0.0340] + [Expression level of CSF3R * 0.1133] + [Expression level of IL2RG * 0.0138] + [Expression level of PTGER4 * 0.1130]. This analysis showed that high risk scores corresponded with high fatality (Figure 4D). Univariate and multivariate Cox regression analyses revealed that the risk score could independently predict the prognosis of patients after adjusting for other parameters like distant metastasis, lymph node metastasis, tumor stage, clinical stage, pathologic grade, age, and gender (Figure 4E, 4F). Thus, the IRGs-based prognostic index may help stratify HCC patients based on expected clinical outcomes.
Figure 4. Clinical utility of prognostic IRG panel and IRGs-based prognostic index. (A) Overall survival, disease-free survival, cancer-specific survival, and progression-free survival of patients between high- and low-risk groups based on prognostic IRGs. (B) ROC curve analysis of risk score compared with clinicopathological features (age, gender, pathologic grade, clinical stage, tumor stage, lymph node metastasis and distant metastasis). (C, D) Differential risk scores, survival status and expression pattern of 6 IRGs in HCC patients. (E, F) Univariate and multivariate Cox regression analysis showing the independent prognostic value of this risk score. (G) Relationships between genes in the prognostic IRG panel and clinicopathological characteristics of HCC patients. (H) Pearson correlation analysis between risk score and infiltration abundances of immune cells. (I) Spearman correlation analysis between risk score and tumor mutation burden (TMB). (J) Spearman correlation analysis between risk score and microsatellite instability (MSI). (K) Heatmap visualization of the expression of 4 DNA mismatch repair (MMR) genes related to risk score. *** p = <0.001.
Clinical utility of the prognostic IRG panel and IRGs-based prognostic index
In further analysis, we examined the relationship between genes in the prognostics panel and clinicopathological features of HCC patients. This analysis revealed a significant negative correlation between IL18RAP and distant metastasis, advanced T stage and advanced clinical stage (Figure 4G, Supplementary Figure 1A). FCER1G exhibited significant positive correlation with poor pathologic grade (Figure 4G). CXCL11 showed significant negative correlation with lymph node metastasis and distant metastasis (Supplementary Figure 1A). CSF3R exhibited a marked negative correlation with distant metastasis (Supplementary Figure 1A). IL2RG exhibited significant negative correlation with lymph node metastasis (Supplementary Figure 1A). To determine if the IRGs-based prognostic index accurately reflects tumor immune microenvironment status, we examined the association of risk scores with infiltration abundance of various types of immune cells. Significant positive correlation was observed between IRGs-based prognostic index and infiltration abundance of memory B-cells or gamma delta T-cells (r=0.527, p=<0.001; r=0.292, p=<0.05) (Figure 4H). Negative correlation was observed between the IRGs-based prognostic index and infiltration abundance of resting memory CD4 T-cells (r=-0.314, P<0.05) (Figure 4H).
Tumor mutation burden (TMB), impaired DNA mismatch repair (MMR), and microsatellite instability (MSI) have been found to affect response to cancer immunotherapy. Interestingly, we find that the prognostic IRG panel-based risk score correlates positively with TMB and MSI (Figure 4I, 4J). Relative to the low risk group, PMS2, MLH1, MSH2 and MSH6 (the MMR genes) were significantly elevated in the high-risk group (Figure 4K). These results strongly suggest that patients in high risk group may benefit from immunotherapy.
Changes in TIME and prognosis based on differential stromal cell infiltration
Re-subtyping HCC according to differential stromal components
Although ssGSEA-based HCC subtyping showed that Immunity_H contains more stromal and immune cells (Figure 1B, 1C), a more rigorous grouping method based on ESTIMATE scores (stromal score combined with immune score) is necessary to determine disease features related to stromal contents. Thus, we reclassified the 369 HCC samples into high stromal content group (Stroma_H) and low stromal content group (Stroma_L) based on median ESTIMATE score. Expectedly, survival analysis suggested that Stroma_H correlated with significantly better 3-year survival than Stroma_L (Figure 5A). To identify the stromal components responsible for distinct clinical outcomes, we further divided the samples into the high stromal cell content group (Stroma cell_H), low stromal cell content group (Stroma cell_L), high immune cell content group (Immune cell_H) and low immune cell content group (Immune cell_L), based on median value of stromal score and immune score, respectively. Survival analyses indicated that Stroma cell_H significantly correlated with better 3-year survival, while the effect of immune cell content on survival rate did not significantly differ between groups (Figure 5B, 5C).
Figure 5. Identification of stromal cell content-related prognostic genes. (A–C) Comparison of 3-year survival between Stroma_H and Stroma_L, Stroma cell_H and Stroma cell_L, and Immune cell_H and Immune cell_L. (D, E) Heatmap and volcano plot of DEGs between Stroma cell_H and Stroma cell_L. (F) Forest plot of hazard ratios and corresponding 95% confidence intervals estimated from univariate Cox regression analyses. Variables significantly associated with good and poor OS are shown in green and red, respectively. (G) Heatmap of the differentially expressed TF genes between Stroma cell_H and Stroma cell_L. (H) Combinatorial TF regulatory networks.
Identification of stromal cell content-related prognostic genes
To identify stromal cell content-related differentially expressed genes (DEGs), we compared Stroma cell_H gene expression profiles to those of Stroma cell_L and identified 480 significantly associated genes. Of these, 455 genes were elevated and 25 genes were suppressed in Stroma cell_H (Figure 5D, 5E). 35 prognostic DEGs significantly associated with favorable or unfavorable HCC survival outcomes were then identified from the above DEGs. A forest plot of hazard ratios showed that 2 of these DEGs were protective factors and 33 were significantly correlated with poor prognosis (Figure 5F). To explore the potential regulatory mechanisms of our prognostic DEGs, we evaluated the expression profiles of the 318 TFs and found 3 to be significantly elevated in Stroma cell_H relative to Stroma cell_L (Figure 5G). The regulatory network analysis showed relationships between 2 of the 3 TFs and 18 of the 35 prognostic DEGs (Figure 5H).
Evaluation of clinical outcomes based on prognostic DEG panel
Next, we assessed the risk score for stromal cell content-related prognostic genes and designed a prognostic signature to subcategorize HCC patients into 2 groups based on prognosis. This analysis found that cancer specific survival, disease-free survival, overall survival, and progression-free survival were higher in the low risk group than in the high risk group (Figure 6A). The AUC of the survival-dependent ROC curve was 0.752, which is higher than the AUCs of clinicopathologic factors, suggesting high prognostic value for this DEG panel (Figure 6B). Heatmaps were used to visualize expression profiles of the 17 genes in the panel from low risk to high risk group (Figure 6C). Risk score was calculated as follows: [Expression level of IL7R * (-0.2979)] + [Expression level of FCGR1A * (-0.2721)] + [Expression level of GAL3ST4 * (-0.2270)] + [Expression level of COL16A1 * (-0.0941)] + [Expression level of ITGA3 * (-0.0522)] + [Expression level of MMP9 * (-0.0144) + [Expression level of CTSV * 0.0166]+[Expression level of MEP1A* 0.0177]+[Expression level of GBP5 * 0.0192]+[Expression level of TMEM130 * 0.0474]+[Expression level of VSIG4 * 0.0873]+[Expression level of CYBB * 0.0882]+[Expression level of HTRA3 * 0.1048]+[Expression level of PLAC8 * 0.1050]+[Expression level of OLFM1 * 0.1162]+[Expression level of FAM183A * 0.1604]+[Expression level of PNMA2 * 0.4091]. High risk scores correlated with higher contemporaneous deaths (Figure 6D). Furthermore, the risk score exhibited independent prognostic potential upon adjustment of other parameters, like distant metastasis, lymph node metastasis, tumor stage, clinical stage, pathologic grade, gender and age (Figure 6E, 6F). Stromal cell content-related, DEG-based prognostic index, exhibited significant positive correlation with advanced clinical stage and advanced T stage (Figure 6G).
Figure 6. Clinical utility of stromal cell content-related prognostic DEG panel and stromal score. (A) Overall survival, disease-free survival, cancer-specific survival, and progression-free survival of patients between high- and low-risk groups based on stromal cell content-associated prognostic genes. (B) ROC curve analysis of risk score relative to clinicopathological features (age, gender, pathologic grade, clinical stage, tumor stage, lymph node metastasis and distant metastasis). (C, D) Differential risk scores, survival status and expression pattern of 17 IRGs in HCC patients. (E, F) Univariate and multivariate Cox regression analysis showing the independent prognostic value of this risk score. (G) Relationship between risk score and the clinical stage, as well as the T stage of HCC patients. (H) Relationships between genes in the prognostic DEG panel and the clinicopathological characteristics of HCC patients. (I) Pearson correlation analysis between the stromal score and infiltration abundances of immune cells. (J) Spearman correlation analysis between the risk score and TMB. (K) Spearman correlation analysis between the risk score and MSI. (L) Heatmap visualization of the expression of 4 DNA MMR genes related to the risk score. *** p = <0.001.
Clinical utility of the prognostic DEG panel and stromal score
Analysis of the relationships between the above genes in DEG panel and HCC clinicopathological features revealed that GBP5, IL7R, PNMA2, ITGA3, COL16A1, FCGR1A, MEP1A, CTSV, FAM183A and PLAC8 expression were strongly linked to distant metastasis, lymph node metastasis, T stage, clinical stage or pathologic grade (Figure 6H, Supplementary Figure 1B). Moreover, analysis of the relationships between stromal cell content-related DEGs-based prognostic index and infiltration abundance of various immune cell types revealed a significant positive correlation between stromal cell content-related, DEG-based prognostic index and infiltration by monocytes, M2 macrophages, or activated mast cells (r=0.353, p=<0.05; r=0.544, p=<0.001; r=0.576, p=<0.001) (Figure 6I). A negative correlation was observed between the stromal cell content-related, DEGs-based prognostic index, and infiltration abundances of regulatory T-cells (r=-0.314, p=<0.05) (Figure 6I). Furthermore, we found that risk score based on the stromal cell content-related prognostic genes positively correlated with MSI (Figure 6K), but did not significantly correlate with TMB (Figure 6J). Relative to the low risk group, expression of the MMR genes PMS2, MLH1, MSH2, and MSH6, was higher in the high risk group (Figure 6L), indicating that patients in the high-risk group derived much benefit from immunotherapy.
Changes in TIME and prognosis based on differential immune cell composition
Identification the heterogeneous of tumor-infiltrating immune cells
Next, we used CIBERSORT [22], to investigate immune infiltration of 22 subpopulations of immune cells in the above HCC samples. Fifty cases with reliably predicted results were identified using a p-value cutoff of ≤0.05. The proportions of immune cells in HCC vary significantly between samples (Figure 7A). Since we speculated that variation in the proportion of tumor-infiltrating immune cells (TIICs) may represent differences between subtypes, the samples of Immunity_H (34/50 cases) and Immunity_L (16/50cases), were separated into 2 discrete groups for further investigation (Figure 7B). Relative to Immunity_L, naive B-cells, memory B-cells, plasma cells, CD8 T-cells, activated CD4 memory T-cells, follicular helper T-cells, gamma delta T-cells, M0 macrophages, M1 macrophages and activated dendritic cells, presented significant abundance differences in Immunity_H (Figure 7C). The proportions of different TIIC subpopulations exhibited only weakly-moderately correlation with each other (Figure 7D). These results indicate the clinical significance of heterogeneous immune infiltration in HCC.
Figure 7. Evaluation of clinical outcomes based on differential TIICs. (A) Relative proportions of 22 immune cell subpopulations in HCC patients. (B) Heatmap visualization of differential immune cell proportions between Immunity_H and Immunity_L. (C) Violin plot analysis exhibiting distinct immune cells subpopulation between Immunity_H and Immunity_L. (D) Correlation matrix of all 22 immune cell proportions in HCC. (E) Box plot depicting relationships between TIICs and pathologic grade, as well as the clinical stage (Kruskal-Wallis tests). (F) Kaplan-Meier survival curves show the relationship between TIICs and survival. (G) The regulation networks among prognostic TIICs, IRG panel, and stromal cell content-related DEG panel. Green and blue represent protective factors. Red and orange represent the opposite.
Evaluation of clinical outcomes based on differential TIICs
The presence of tumor-infiltrating lymphocytes (TILs) in TIICs is associated with better responses to immunotherapy and improved clinical outcomes [23]. Expectedly, resting dendritic cells were strongly associated with higher pathologic grade and lower stage, whereas M1 macrophages and CD8 T-cells strongly related to lower pathologic grade (Figure 7E). Survival analyses revealed that higher proportions of TILs (resting memory CD4 T-cells and activated NK cells), and lower proportions of eosinophils and activated mast cells, significantly correlated with better 3- or 5-year survival rates (Figure 7F). A TIICs-based regulatory schematic was used to illustrate potential interactions between prognostic TIICs, IRG panel, and stromal cell content-related DEG panel. Most poor prognosis-associated genes (12/18) correlated negatively with favorable prognostic TILs (resting memory CD4 T-cells and activated NK cells) (Figure 7G). We speculate these prognostic genes remodel TIME, leading to differential immune cell infiltration, thereby influencing HCC clinical outcomes.
Discussion
Multiple studies have subtyped HCC using genomics, transcriptomics, and metabolomics [24–26]. However, few have examined HCC classification by immune signatures. Here, we identified immune-related HCC subtypes using immunogenomic profiles, stromal cell features, and immune cell composition. Our data show that intratumoral genetic and immune microenvironment heterogeneities are essential features of HCC. Understanding individualized immune signature as a potential avenue is expected to improve immunotherapeutic responsiveness.
Tumor gene expression profiling has identified gene expression signatures with prognostic value that can inform patient stratification for targeted therapies. Recent studies have evaluated the expression pattern of IRGs in solid tumors patients receiving immune-based therapies. For instance, in samples from melanoma patients receiving recombinant IL2 treatment, a signature that could predict the clinical response was identified [27]. Elsewhere, an IFN-inflammatory immune gene expression-based signature, which correlated with progression-free survival and enhanced overall response rates, was established in melanoma patients receiving pembrolizumab. The signature is also being investigated in other malignancies [28]. T-effector/IFNγ signature, a 8-gene signature reflecting preexisting immunity, is in phase II trial on previously treated non-small cell lung carcinoma (NSCLC) [29]. Here, we classified HCC into 2 stable subtypes (Immunity High and Immunity Low) based on ssGSEA score. From 729 genes exhibiting differential expression between subtypes, a panel of 6 IRGs was found to significantly correlate with survival outcome, with high prognostic value. Moreover, expression of individual genes in this panel, including IL18RAP, CXCL11, FCER1G, CSF3R and IL2RG were strongly linked to distant metastasis, lymph node metastasis, tumor stage, clinical stage or pathologic grade. In Immunity_H subtype, HLA genes, PD-L1 protein, chemotactic factors related to immunologic activity, MHC-class II protein complex, natural killer cell chemotaxis, antigen processing and presentation were upregulated. This could lead to favorable clinical outcomes. The positive correlation between the prognostic IRG panel-based risk score and TMB, MSI or MMR genes expression also strongly suggested that patients in the high risk group might benefit more from immunotherapy. Thus, identification of intratumoral heterogeneity with regards to prognostic IRG panel, immune-related cellular component, biological process and pathways are used to build prognostic models of pathologic grade, TNM stages, and survival. It may also guide the selection of potential immunotherapeutic targets and individualize treatments.
Tumors possess a stromal compartment comprising cellular and non-cellular components, which also influence cancer development, progression, and metastasis [30, 31]. Although studies have shown that anti-cancer therapies target cancer cells, their effect on the tumor stroma has not been well defined. As the 2 main types of non-tumor components in tumor immune microenvironment, stromal cells and immune cells are considered necessary in cancer diagnostic and prognostic applications [32]. In bulk urothelial cancer transcriptomes, Li Wang et al. recently reported that the primary source of EMT-related gene expression is the non-hematopoietic stromal cells. In a study sample comprising patients with metastatic urothelial cancer under nivolumab treatment, higher EMT/stroma-associated genes expression were related to lower response rate, shorter progression-free survival and poorer overall survival in patient tumors with T-cell infiltration [33]. Here, we grouped HCC into high stromal cell content group (Stroma cell_H), low stromal cell content group (Stroma cell_L), high immune cell content group (Immunecell_H) and low immune cell content group (Immune cell_L) based on stromal score or immune score and evaluated stromal component-related pathogenesis. Our data show that Stroma cell_H correlates with improved 3-year survival significantly. A panel of 17 stromal cell content related DEGs between Stroma cell_H and Stroma cell_L was found to correlate with survival rate significantly. Moreover, the expression level of individual genes in this panel, including GBP5, IL7R, PNMA2, ITGA3, COL16A1, FCGR1A, MEP1A, CTSV, FAM183A and PLAC8, were strongly associated with distant metastasis, lymph node metastasis, tumor stage, clinical stage, or pathologic grade. A positive correlation between risk score based on stromal cell content-related prognostic genes and MSI or MMR genes expression also strongly showed that the patients in the high risk group are more likely to benefit from immunotherapy. These data reflect high heterogeneity of tumor stromal cells in HCC, suggesting that anti-cancer therapies should not only target cancer cells but also the stromal compartment for effective outcomes.
The terms ‘cold’ and ‘hot’ refer to non-inflamed tumors, inflamed but non-infiltrated, and T cell-infiltrated, reflecting lower and higher Immunoscore categories [34]. Hot tumors are characterized by the presence of TILs, expression of anti-PD-L1 on tumor-associated immune cells, potential genomic instabilities and existence of a previous anti-tumor immune response [35–37]. These tumors had a higher response to immune-based therapy. Here, we show that multiple immune cell types were present at significantly higher abundance in Immunity_H, especially the TILs, including memory B-cells, naive B-cells, activated CD4 memory T-cells, CD8 T-cells, follicular helper T cells, and gamma delta T-cells. Higher proportions of TILs significantly correlated with improved 3-, 5- or 10-year survival rates. Moreover, Immunity_H displayed significantly elevated PD-L1 level relative to Immunity_L. Given that TILs are relatively rare in HCC [38], our study highlights immune cell heterogeneity in HCC, suggesting a novel way of stratifying patients for immunotherapy. Our TIICs-based regulatory schematic also highlights avenues for converting immune cold tumors into hot tumors.
Immunophenotypes of solid tumors fall into 3 subtypes: inflamed, immune excluded, and immune desert [32]. 3 distinct immune-cell infiltration subtypes have recently been identified [39]. Relative to previously reported immune subtypes, our clustering based on ssGSEA analysis identified only 2 subtypes. This may arise from the fact that our study focus on the prognostic value. The two subtypes, Immunity_H group and Immunity_L group, exhibited distinct survival rates in both the TCGA and GSE14520 datasets (Figure 1H). However, the prognostic value of Immunity_M group was not consistent across datasets when samples were grouped into 3 groups (data not shown), indicating that 2 subtypes are more effective and that Immunity_H group predicts better survival. In some studies, classification was mainly based on immune-cell infiltration. Here, we took both immune-cell infiltration and immune-related pathways into consideration. This is because high immune-cell infiltration does not necessarily mean high immune response because cells may not be activated. Thus, our strategy may more accurately reflect the immune microenvironment landscape.
Given the significance of HCC immune landscape, our study improves our understanding of intratumoral genetic and immune microenvironment heterogeneity from tumor and non-tumor components perspectives. Here, we reveal potential molecular mechanisms and approaches to manipulate the immune status to improve individualized immunotherapy.
Materials and Methods
Data acquisition and clinical samples
FPKM-normalized RNA-sequencing data from 369 primary HCC cases and corresponding prognostic data were retrieved from TCGA (https://https://portal.gdc.cancer.gov).
Clustering
Enrichment level and activity of several immune cells, pathways or functions in HCC were analyzed using single sample gene-set enrichment analysis (ssGSEA) score based on 29 immune-associated gene sets [12]. Hierarchical clustering of HCC was done based on the ssGSEA scores of the 29 immune signatures using the “sparcl” package on R.
Estimation of immune score, stromal score, and tumor purity
The normalized expression matrix was analyzed by ESTIMATE and the immune score (the infiltration level of immune cells), stromal score (the level of stromal cells present) and tumor purity for each HCC sample determined [40].
Identification of differentially expressed genes
Differentially expressed genes were identified using “limma”, “ggpubr” and “pheatmap” packages, with FDR =<0.05 and a log2 |fold change| >1 as cutoffs.
Gene-set enrichment analysis
Functional enrichment analyses via GO term and KEGG pathway analyses, were performed using GSEA (R implementation) [41–43].
Identification of immune-related gene
IRGs were identified using ImmPort (immunology database and analysis portal) database [44]. These genes are known to be involved in immune activity. Differentially expressed IRGs were extracted from all differentially expressed genes.
Identification of survival-associated genes and survival analyses
Survival-associated genes were selected using univariate Cox’s proportional hazards regression analysis. Kaplan-Meier analysis was used to compare survival differences. Log-rank test was used to calculate statistical significance with p = <0.05 as the threshold.
Identification of gene-based prognostic index (GPI)
Survival-associated genes were subjected to multivariate Cox’s regression analyses, with integrated genes panels as independent prognostic indicators for GPI development. Immune-related gene-based prognostic index (IRGPI) was based on expression data multiplied by the Cox regression coefficient. Patients were divided into high- and low-risk groups based on the median risk score.
Construction of transcription factors (TFs) regulatory network
A list of 318 TFs was obtained from the Cistrome cancer database which integrates cancer genomics data from TCGA with over 23000 ChIP-seq chromatin profiling data from Cistrome, illuminating regulatory links between TFs and the transcriptome [45].
Evaluation of the tumor-infiltrating immune cells (TIICs) proportion
Normalized gene expression data were used to estimate the fraction of 22 infiltrating immune cell types using CIBERSORT as previously described [22]. 1000 CIBERSORT permutations and cases were set with p= <0.05 as a cutoff value. Mann-Whitney U test was used to compare the proportion of immune cell subsets between HCC subtypes.
TMB analysis
Somatic mutation data on 357 hepatocellular carcinoma patients were obtained and processed using VarScan software from the “Masked Somatic Mutation” category on TCGA. Next, TMB was defined and calculated using the formula: (total count of variants)/(the whole length of exons). Detected variants included base substitutions, deletions, or insertions.
MSI analysis
MSI status for 367 HCC cases consecutively sequenced with Memorial Sloan Kettering-Integrated Mutation Profiling of Actionable Cancer Targets clinical NGS assay was determined using MSI sensor algorithm, a program that reports the percentage of unstable microsatellites [46].
Statistical analysis
To examine performance of the prognostic index, time-dependent receiver operating characteristic (ROC) curve was calculated using “survivalROC” package [47]. Multivariate Cox regression analyses were performed to verify if the risk score was an independent prognostic index. All analyses were performed on R and the regulatory network visualized using Cytoscape software version 3.7.1 (https://cytoscape.org/). Statistical tests were two-sided. p = <0.05 indicated statistical significance.
Supplementary Materials
Supplementary Figure 1
Author Contributions
Hongmei Yan and Yiyi Li: Study design, implementation and drafting. Yuchuan Chen and Kai Wang: Data acquisition and analysis. Lu Yu and Yuwen Xie: Data collection and organization. Xixin Huang, Qianyu Li, Jiayu Lin, Yueyun He, Xinyu Yi and Yanzhi Wang: Data interpretation. Xixin Huang, Longhua Chen and Yi Ding: Manuscript revise and final review.
Conflicts of Interest
The authors have no conflicts of interests to declare.
Funding
This work is supported by National Natural Science Foundation of China (No. 81803063, No. 81872470 and No. 81672992), Natural Science Foundation of Guangdong Province (No. 2018030310297), Guangdong Natural Science Funds for Distinguished Young Scholar (No. 2015A030306015), Guangdong Program for Support of Top-notch Young Professionals (No. 2015TQ01R279), Outstanding Youths Development Scheme of Nanfang Hospital, Southern Medical University (No. 2018J009), Foundation of President of Nanfang Hospital, Southern Medical University (No. 2017C016), and Provincial College Students' Innovation and Entrepreneurship training Program, Southern Medical University (No. S201912121083 and No. S201912121079).
References
- 1. Crissien AM, Frenette C. Current management of hepatocellular carcinoma. Gastroenterol Hepatol (N Y). 2014; 10:153–61. [PubMed]
- 2. Buonaguro L, Mauriello A, Cavalluzzo B, Petrizzo A, Tagliamonte M. Immunotherapy in hepatocellular carcinoma. Ann Hepatol. 2019; 18:291–97. https://doi.org/10.1016/j.aohep.2019.04.003 [PubMed]
- 3. Zhu XD, Sun HC. Emerging agents and regimens for hepatocellular carcinoma. J Hematol Oncol. 2019; 12:110. https://doi.org/10.1186/s13045-019-0794-6 [PubMed]
- 4. Fu Y, Liu S, Zeng S, Shen H. From bench to bed: the tumor immune microenvironment and current immunotherapeutic strategies for hepatocellular carcinoma. J Exp Clin Cancer Res. 2019; 38:396. https://doi.org/10.1186/s13046-019-1396-4 [PubMed]
- 5. Xu W, Liu K, Chen M, Sun JY, McCaughan GW, Lu XJ, Ji J. Immunotherapy for hepatocellular carcinoma: recent advances and future perspectives. Ther Adv Med Oncol. 2019; 11:1758835919862692. https://doi.org/10.1177/1758835919862692 [PubMed]
-
6.
El-Khoueiry AB, Sangro B, Yau T, Crocenzi TS, Kudo M, Hsu C, Kim TY, Choo SP, Trojan J, Welling TH
Rd , Meyer T, Kang YK, Yeo W, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet. 2017; 389:2492–502. https://doi.org/10.1016/S0140-6736(17)31046-2 [PubMed] - 7. Hinohara K, Polyak K. Intratumoral heterogeneity: more than just mutations. Trends Cell Biol. 2019; 29:569–79. https://doi.org/10.1016/j.tcb.2019.03.003 [PubMed]
- 8. Li L, Li M, Jiang Z, Wang X. ARID1A mutations are associated with increased immune activity in gastrointestinal cancer. Cells. 2019; 8:678. https://doi.org/10.3390/cells8070678 [PubMed]
- 9. Dong ZY, Zhang C, Li YF, Su J, Xie Z, Liu SY, Yan LX, Chen ZH, Yang XN, Lin JT, Tu HY, Yang JJ, Zhou Q, et al. Genetic and immune profiles of solid predominant lung adenocarcinoma reveal potential immunotherapeutic strategies. J Thorac Oncol. 2018; 13:85–96. https://doi.org/10.1016/j.jtho.2017.10.020 [PubMed]
- 10. Runa F, Hamalian S, Meade K, Shisgal P, Gray PC, Kelber JA. Tumor microenvironment heterogeneity: challenges and opportunities. Curr Mol Biol Rep. 2017; 3:218–29. https://doi.org/10.1007/s40610-017-0073-7 [PubMed]
- 11. Catalano I, Trusolino L. The stromal and immune landscape of colorectal cancer progression during anti-EGFR therapy. Cancer Cell. 2019; 36:1–3. https://doi.org/10.1016/j.ccell.2019.06.001 [PubMed]
- 12. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
- 13. He Y, Jiang Z, Chen C, Wang X. Classification of triple-negative breast cancers based on immunogenomic profiling. J Exp Clin Cancer Res. 2018; 37:327. https://doi.org/10.1186/s13046-018-1002-1 [PubMed]
- 14. Jiang X, Wang J, Deng X, Xiong F, Ge J, Xiang B, Wu X, Ma J, Zhou M, Li X, Li Y, Li G, Xiong W, et al. Role of the tumor microenvironment in PD-L1/PD-1-mediated tumor immune escape. Mol Cancer. 2019; 18:10. https://doi.org/10.1186/s12943-018-0928-4 [PubMed]
- 15. Yarchoan M, Albacker LA, Hopkins AC, Montesion M, Murugesan K, Vithayathil TT, Zaidi N, Azad NS, Laheru DA, Frampton GM, Jaffee EM. PD-L1 expression and tumor mutational burden are independent biomarkers in most cancers. JCI Insight. 2019; 4:e126908. https://doi.org/10.1172/jci.insight.126908 [PubMed]
- 16. Ruiz de Galarreta M, Bresnahan E, Molina-Sánchez P, Lindblad KE, Maier B, Sia D, Puigvehi M, Miguela V, Casanova-Acebes M, Dhainaut M, Villacorta-Martin C, Singhi AD, Moghe A, et al. Β-catenin activation promotes immune escape and resistance to anti-PD-1 therapy in hepatocellular carcinoma. Cancer Discov. 2019; 9:1124–41. https://doi.org/10.1158/2159-8290.CD-19-0074 [PubMed]
- 17. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010; 26:139–40. https://doi.org/10.1093/bioinformatics/btp616 [PubMed]
- 18. Kreiter S, Vormehr M, van de Roemer N, Diken M, Löwer M, Diekmann J, Boegel S, Schrörs B, Vascotto F, Castle JC, Tadmor AD, Schoenberger SP, Huber C, et al. Mutant MHC class II epitopes drive therapeutic immune responses to cancer. Nature. 2015; 520:692–96. https://doi.org/10.1038/nature14426 [PubMed]
- 19. Jin Y, Knudsen E, Wang L, Maghazachi AA. Lysophosphatidic acid induces human natural killer cell chemotaxis and intracellular calcium mobilization. Eur J Immunol. 2003; 33:2083–89. https://doi.org/10.1002/eji.200323711 [PubMed]
- 20. Kondo K, Takada K, Takahama Y. Antigen processing and presentation in the thymus: implications for T cell repertoire selection. Curr Opin Immunol. 2017; 46:53–57. https://doi.org/10.1016/j.coi.2017.03.014 [PubMed]
- 21. Jiang GM, Wang HS, Du J, Ma WF, Wang H, Qiu Y, Zhang QG, Xu W, Liu HF, Liang JP. Bortezomib relieves immune tolerance in nasopharyngeal carcinoma via STAT1 suppression and indoleamine 2,3-dioxygenase downregulation. Cancer Immunol Res. 2017; 5:42–51. https://doi.org/10.1158/2326-6066.CIR-16-0102 [PubMed]
- 22. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–57. https://doi.org/10.1038/nmeth.3337 [PubMed]
- 23. Beyrend G, van der Gracht E, Yilmaz A, van Duikeren S, Camps M, Höllt T, Vilanova A, van Unen V, Koning F, de Miranda NF, Arens R, Ossendorp F. PD-L1 blockade engages tumor-infiltrating lymphocytes to co-express targetable activating and inhibitory receptors. J Immunother Cancer. 2019; 7:217. https://doi.org/10.1186/s40425-019-0700-3 [PubMed]
- 24. Cancer Genome Atlas Research Network. Electronic address: wheeler@bcm.edu, and Cancer Genome Atlas Research Network. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell. 2017; 169:1327–41.e23. https://doi.org/10.1016/j.cell.2017.05.046 [PubMed]
- 25. Kwee SA, Tiirikainen M, Sato MM, Acoba JD, Wei R, Jia W, Le Marchand L, Wong LL. Transcriptomics associates molecular features with 18F-fluorocholine PET/CT imaging phenotype and its potential relationship to survival in hepatocellular carcinoma. Cancer Res. 2019; 79:1696–704. https://doi.org/10.1158/0008-5472.CAN-18-3837 [PubMed]
- 26. Jiang Y, Sun A, Zhao Y, Ying W, Sun H, Yang X, Xing B, Sun W, Ren L, Hu B, Li C, Zhang L, Qin G, et al, and Chinese Human Proteome Project (CNHPP) Consortium. Proteomics identifies new therapeutic targets of early-stage hepatocellular carcinoma. Nature. 2019; 567:257–61. https://doi.org/10.1038/s41586-019-0987-8 [PubMed]
-
27.
Weiss GR, Grosh WW, Chianese-Bullock KA, Zhao Y, Liu H, Slingluff CL
Jr , Marincola FM, Wang E. Molecular insights on the peripheral and intratumoral effects of systemic high-dose rIL-2 (aldesleukin) administration for the treatment of metastatic melanoma. Clin Cancer Res. 2011; 17:7440–50. https://doi.org/10.1158/1078-0432.CCR-11-1650 [PubMed] - 28. Ribas A, Robert C, Hodi FS, Wolchok JD, Joshua AM, Hwu W, Weber JS, Zarour HM, Kefford R, Loboda A, Albright A, Kang SP, Ebbinghaus S, et al. Association of response to programmed death receptor 1 (PD-1) blockade with pembrolizumab (MK-3475) with an interferon-inflammatory immune gene signature. J Clin Oncol. 2015. https://doi.org/10.1200/jco.2015.33.15_suppl.3001
- 29. Fehrenbacher L, Spira A, Ballinger M, Kowanetz M, Vansteenkiste J, Mazieres J, Park K, Smith D, Artal-Cortes A, Lewanski C, Braiteh F, Waterkamp D, He P, et al, and POPLAR Study Group. Atezolizumab versus docetaxel for patients with previously treated non-small-cell lung cancer (POPLAR): a multicentre, open-label, phase 2 randomised controlled trial. Lancet. 2016; 387:1837–46. https://doi.org/10.1016/S0140-6736(16)00587-0 [PubMed]
- 30. Wen Q, Xiong W, He J, Zhang S, Du X, Liu S, Wang J, Zhou M, Ma L. Fusion cytokine IL-2-GMCSF enhances anticancer immune responses through promoting cell-cell interactions. J Transl Med. 2016; 14:41. https://doi.org/10.1186/s12967-016-0799-7 [PubMed]
- 31. Zhang M, Yan Z, Wang J, Yao X. Toll-like receptors 7 and 8 expression correlates with the expression of immune biomarkers and positively predicts the clinical outcome of patients with melanoma. Onco Targets Ther. 2017; 10:4339–46. https://doi.org/10.2147/OTT.S136194 [PubMed]
- 32. Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, Coussens LM, Gabrilovich DI, Ostrand-Rosenberg S, Hedrick CC, Vonderheide RH, Pittet MJ, Jain RK, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med. 2018; 24:541–50. https://doi.org/10.1038/s41591-018-0014-x [PubMed]
- 33. Wang L, Saci A, Szabo PM, Chasalow SD, Castillo-Martin M, Domingo-Domenech J, Siefker-Radtke A, Sharma P, Sfakianos JP, Gong Y, Dominguez-Andres A, Oh WK, Mulholland D, et al. EMT- and stroma-related gene expression and resistance to PD-1 blockade in urothelial cancer. Nat Commun. 2018; 9:3503. https://doi.org/10.1038/s41467-018-05992-x [PubMed]
- 34. Sweis RF, Spranger S, Bao R, Paner GP, Stadler WM, Steinberg G, Gajewski TF. Molecular drivers of the non-T-cell-inflamed tumor microenvironment in urothelial bladder cancer. Cancer Immunol Res. 2016; 4:563–68. https://doi.org/10.1158/2326-6066.CIR-15-0274 [PubMed]
- 35. Teng MW, Ngiow SF, Ribas A, Smyth MJ. Classifying cancers based on T-cell infiltration and PD-L1. Cancer Res. 2015; 75:2139–45. https://doi.org/10.1158/0008-5472.CAN-15-0255 [PubMed]
- 36. Mouw KW, Goldberg MS, Konstantinopoulos PA, D’Andrea AD. DNA damage and repair biomarkers of immunotherapy response. Cancer Discov. 2017; 7:675–93. https://doi.org/10.1158/2159-8290.CD-17-0226 [PubMed]
- 37. Sun XG, Lin XC, Diao JX, Yu ZL, Li K. Pi (spleen)-deficiency syndrome in tumor microenvironment is the pivotal pathogenesis of colorectal cancer immune escape. Chin J Integr Med. 2016; 22:789–94. https://doi.org/10.1007/s11655-015-2086-5 [PubMed]
- 38. Shirabe K, Motomura T, Muto J, Toshima T, Matono R, Mano Y, Takeishi K, Ijichi H, Harada N, Uchiyama H, Yoshizumi T, Taketomi A, Maehara Y. Tumor-infiltrating lymphocytes and hepatocellular carcinoma: pathology and clinical management. Int J Clin Oncol. 2010; 15:552–8. https://doi.org/10.1007/s10147-010-0131-0 [PubMed]
- 39. Kurebayashi Y, Ojima H, Tsujikawa H, Kubota N, Maehara J, Abe Y, Kitago M, Shinoda M, Kitagawa Y, Sakamoto M. Landscape of immune microenvironment in hepatocellular carcinoma and its additional impact on histological and molecular classification. Hepatology. 2018; 68:1025–41. https://doi.org/10.1002/hep.29904 [PubMed]
- 40. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
- 41. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017; 45:D353–61. https://doi.org/10.1093/nar/gkw1092 [PubMed]
- 42. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102 [PubMed]
- 43. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
- 44. Bhattacharya S, Andorf S, Gomes L, Dunn P, Schaefer H, Pontius J, Berger P, Desborough V, Smith T, Campbell J, Thomson E, Monteiro R, Guimaraes P, et al. ImmPort: disseminating data to the public for the future of immunology. Immunol Res. 2014; 58:234–39. https://doi.org/10.1007/s12026-014-8516-1 [PubMed]
- 45. Mei S, Meyer CA, Zheng R, Qin Q, Wu Q, Jiang P, Li B, Shi X, Wang B, Fan J, Shih C, Brown M, Zang C, Liu XS. Cistrome cancer: a web resource for integrative gene regulation modeling in cancer. Cancer Res. 2017; 77:e19–22. https://doi.org/10.1158/0008-5472.CAN-17-0327 [PubMed]
- 46. Middha S, Zhang L, Nafa K, Jayakumaran G, Wong D, Kim HR, Sadowska J, Berger MF, Delair DF, Shia J, Stadler Z, Klimstra DS, Ladanyi M, et al. Reliable pan-cancer microsatellite instability assessment by using targeted next-generation sequencing data. JCO Precis Oncol. 2017; 2017:10.1200. https://doi.org/10.1200/PO.17.00084 [PubMed]
- 47. Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000; 56:337–44. https://doi.org/10.1111/j.0006-341x.2000.00337.x [PubMed]