Shanshan Liu
Guangchuang Yu
Li Liu
Xuejing Zou1,2- 1Country Guangdong Provincial Key Laboratory of Viral Hepatitis Research, Hepatology Unit and Department of Infectious Diseases, Nanfang Hospital, Southern Medical University, Guangzhou, China
- 2State Key Laboratory of Organ Failure Research, Nanfang Hospital, Southern Medical University, Guangzhou, China
- 3Department of Bioinformatics, School of Basic Medical Sciences, Southern Medical University, Guangzhou, China
- 4Department of Medical Quality Management, Nanfang Hospital, Southern Medical University, Guangzhou, China
- 5Department of Radiation Oncology, Nanfang Hospital, Southern Medical University, Guangzhou, China
A growing amount of evidence has suggested the clinical importance of stromal and immune cells in the liver cancer microenvironment. However, reliable prognostic signatures based on assessments of stromal and immune components have not been well-established. This study aimed to identify stromal-immune score–based potential prognostic biomarkers for hepatocellular carcinoma. Stromal and immune scores were estimated from transcriptomic profiles of a liver cancer cohort from The Cancer Genome Atlas using the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumors using Expression data) algorithm. Least absolute shrinkage and selection operator (LASSO) algorithm was applied to select prognostic genes. Favorable overall survivals and progression-free interval were found in patients with high stromal score and immune score, and 828 differentially expressed genes were identified. Functional enrichment analysis and protein–protein interaction networks further showed that these genes mainly participated in immune response, extracellular matrix, and cell adhesion. MMP9 (matrix metallopeptidase 9) was identified as a prognostic tumor microenvironment–associated gene by using LASSO and TIMER (Tumor IMmune Estimation Resource) algorithms and was found to be positively correlated with immunosuppressive molecules and drug response.
Introduction
Hepatocellular carcinoma (HCC) is the third leading cause of cancer death worldwide. The median survival of HCC patients in China is about 23 months, and ≥ 60% of patients present with intermediate-stage or advanced-stage HCC (Kanwal and Singal, 2019; Yang et al., 2019). Currently, the main treatment for HCC patients in early stages is surgery, combination with transarterial chemoembolization, ablation, and liver transplantation. For others in advanced stages, the effective approaches involve molecular targeting agents (sorafenib, lenvatinib, and regorafenib). Although these methods have improved the prognosis of HCC patients, the overall survival (OS) of HCC remains challenging for the heterogeneity of HCC. And also, there is still a lack of molecular markers used in determination of prognosis and treatment for patients (Bruix et al., 2016).
The liver cancer microenvironment consists of not only tumor cells but also stromal cells, including distinct immune cell subsets. Tumor-infiltrating immune cells and stromal cells are associated with angiogenesis, immune suppression, chemotherapeutic resistance, and tumor cell migration (Affo et al., 2017; Barry et al., 2020; Jin and Jin, 2020; Son et al., 2020; Zhang et al., 2020). An increasing amount of evidence has suggested the clinical importance of stromal cells and immune cells in the microenvironment of liver cancer tissues, tumor microenvironment (TME)–associated genes also have potential as novel biomarkers for a range of cancers (Yang et al., 2020).
In the present study, the Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) algorithm (Yoshihara et al., 2013) was applied to estimate the stromal and immune scores of a series of cancer tissues based on their transcriptional profiles, to perform a comprehensive analysis of immune and stromal cells, and to correlate the data to clinical outcomes of patients.
The least absolute shrinkage and selection operator (LASSO) method is a compressed estimation used to obtain a refined model by constructing a penalty function (Korenberg, 2006). It can help with the selection of variables at the time of parameter estimation so as to better solve the multicollinearity problem of regression analysis. A growing body of research confirms that LASSO is an effective method for gene selection of tumors (Wang et al., 2020; Xu et al., 2020).
Tumor IMmune Estimation Resource (TIMER) integrates multiple state-of-the-art algorithms for immune infiltration estimation, which can explore various associations between immune infiltrates and genetic features in The Cancer Genome Atlas (TCGA) cohorts (Li et al., 2017, 2020). Computational Analysis of REsistance (CARE) is a computational method focused on targeted therapies, to infer genome-wide transcriptomic signatures of drug efficacy from cell line compound screens (Jiang et al., 2018). Previous studies have confirmed that the efficacy of immunotherapy is strongly influenced by the composition and abundance of immune cells in the TME (Boyero et al., 2020).
Thus, we combined LASSO, TIMER algorithms, and CARE to preliminarily demonstrate that the expression of TME-associated genes could be new prognostic and reliable drug response biomarkers for HCC patients.
Materials and Methods
Database
In total, data from 365 HCC patients and 18,161 RNAs extracted from RNA-seq data according to ENSEMBL Genomes (hg38) were analyzed in this study. All RNA expression data and the corresponding clinical data were obtained from TCGA (data version, July 19, 2019)1. The clinicopathological characteristics of the analyzed patients are listed in Supplementary Table 1. The progression-free interval (PFI) is characterized as a time without a new tumor occurrence or a death from cancer. The Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) algorithm was applied to the normalized expression matrix for estimating the stromal and immune scores by using “estimate” R package in R software (version: 3.6.3) for each HCC sample.
Correlations Between Prognoses and Stromal/Immune Scores
OS and PFI was used as the primary prognosis endpoint and was estimated by the GraphPad Prism 8.0. Supplementary Figure 3B is realized by R package “Survival” (Therneau, 2020), “Survminer” (Kassambara et al., 2019), and “timeROC” (Paul Blanche and Jacqmin-Gadda, 2013). Based on the stromal and immune scores estimated from each sample, patients were classified into two groups by using X-tile, and prognoses for each group were examined. The bioinformatics tool, X-tile (Camp et al., 2004), was used to determine the optimum cutoff point according to the minimum P-value defined by the Kaplan–Meier analysis and log-rank test. The principle of X-tile is “enumeration method that different values are grouped as truncation values to conduct statistical tests, and the test result with the lowest P-value can be considered as the best truncation value. The survival outcomes of the two groups were compared by log-rank tests. P < 0.05 was considered as statistically significant.
Identification of Differentially Expressed Genes
Data analysis was performed using an open-source web tool NetworkAnalyst2 (Xia et al., 2013a,b; Zhou et al., 2019). Log2 fold change > 1 and adjusted P < 0.05 were set as the cutoffs to screen for differentially expressed genes (DEGs). A website Venn diagrams tool (Bardou et al., 2014)3 was used to identify the commonly upregulated or downregulated DEGs in the immune and stromal groups. Heatmaps and clustering were generated using the R package “ggplot2” (Wickham, 2016), “ggtree” (Yu, 2020b), and “aplot” (Yu, 2020a).
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analyses
GO (Gene Ontology) enrichment analyses were performed by the “Goseq” (Young et al., 2010) R package, and visualization of bubble diagrams used Hiplot4. KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses and visualization of intersection genes were performed by the “clusterProfiler” (Yu et al., 2012) R package and “enrichplot” (Yu, 2019) R package with P < 0.05 as the cutoff value.
Protein–Protein Interaction Network Construction
The protein–protein interaction (PPI) network was retrieved from Search Tool for the Retrieval of Interaction Gene/Proteins (STRING) (Szklarczyk et al., 2019) database with high confidence (0.7) and reconstructed via the Cytoscape software (Shannon et al., 2003). In Cytoscape, we used Molecular COmplex DEtection (MCODE) (Bader and Hogue, 2003) to select two clusters that contained the largest number of nodes. ClueGo (Bindea et al., 2009) App was used to perform enrichment analysis of each cluster selected by MCODE.
Identification of TME-Associated Prognostic Genes
LASSO algorithm was used to identify candidate genes by “glmnet” (Friedman et al., 2010) R package with the number of lambda = 1,000. Clinical outcomes and gene expression profiles were analyzed by LASSO. Lambda.min is the cutoff point that brings minimum mean cross-validated error. Genes with the highest lambda values were selected for further analysis.
Identification of TME-Associated Prognostic Genes
The TIMER algorithm was used to calculate the tumor abundance of six infiltrating immune cells (CD4+ T cells, CD8+ T cells, B cells, neutrophils, macrophages, and dendritic cells) based on RNA-Seq expression profiles data. The correlation between the selected prognostic genes and immune cells was calculated by Spearman correlation analysis by TIMER. The estimation results were calculated by TIMER2.0, CIBERSORT, quanTIseq, xCell, MCP-counter, and EPIC methods. Relations between immunoinhibitors and expression of matrix metallopeptidase 9 (MMP9) were calculated by Spearman correlation analysis by a web tool TISIDB5 (Ru et al., 2019). The correlation coefficient value <0.3 indicates that the correlation is negligible, whereas the correlation coefficient ≥ 0.3 indicates a positive/negative correlation. The CARE software6 was used to identify genome-scale biomarkers of targeted therapy response using compound screen data. For each gene, the CARE score indicates the association between its molecular alteration and drug efficacy. A positive score indicates a higher expression value (or presence of mutation) to be associated with drug response, whereas a negative score indicates drug resistance.
Statistical Analysis
Unpaired t-test was used to compare two groups of continuously distributed variables. Jonckheere–Terpstra test was used to compare three or more groups of continuously distributed variables. The FDR correction was performed in multiple tests. ∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001.
Results
Association of Stromal and Immune Scores With HCC Pathology and Prognosis
A cohort containing 365 liver hepatocellular carcinoma patients with available expression data and clinical information in TCGA database was analyzed. The general pipeline of the data analysis protocol is shown in Figure 1, and the links of tools are listed in Supplementary Table 6. The clinicopathological characteristics of the analyzed patients are listed in Supplementary Table 1. Based on the gene expression data, immune and stromal scores were calculated using the ESTIMATE algorithm (Supplementary Table 2). The associations of stromal and immune scores with HCC patient pathological characteristics were examined by comparing the score distributions among different tumor stages and differentiation grades.
Figure 1. The general pipeline of the data analysis protocol.
Significant associations were observed between stromal scores and tumor differentiation grades; tumors with poorer differentiation yielded higher stromal scores than those differentiated well (Jonckheere–Terpstra test, P = 0.002) (Figures 2A,B).
Figure 2. Relationship between immune and stromal scores and HCC clinical and pathological data. (A,B) Distribution of immune and stromal scores of HCC grades. (C,D) Distribution of immune and stromal scores of AFP value of HCC. AFP is divided into high and low groups at the limit of 400 ng/mL. (E,F) Distribution of immune and stromal scores of new tumor event after initial treatment of HCC. Unpaired t-test was used to compare two groups of continuously distributed variables. Jonckheere–Terpstra test was used to compare three or more groups of continuously distributed variables. ∗P < 0.05 and ∗∗P < 0.01.
As previously described, serum α-fetoprotein (AFP) values are not only of diagnostic value but also of prognostic significance in patients with HCC (Galle et al., 2019). Thus, we compared changes in immune and stromal scores between AFP low (AFP ≤ 400 ng/mL) and high (AFP > 400 ng/mL) samples. The AFP high cases had the lowest stromal scores (unpaired t-test, P = 0.0204) (Figure 2D). Evidence suggests that AFP plays an immune-suppressing role (Yang et al., 2018), but we found that there is no significant difference in the immune score as shown in Figure 2C. We further used the TIMER algorithm to evaluate the effect of AFP on the immune infiltration of HCC, and results showed that the expression of AFP was weakly correlated with the infiltration abundance of the six immune cells (Supplementary Figure 1A). AFP is dynamic in the occurrence and development of HCC, whereas TCGA patients were only tested for AFP at the time of initial diagnosis, which may lead to the bias of the results in our study.
Also, when we compared the immune and stromal scores between patients with a new tumor event and without new tumor event after initial treatment, patients without a new tumor event had higher immune and stromal scores (unpaired t-test, P = 0.0461 for stromal score and p = 0.1966 for immune score) (Figures 2E,F).
We also analyzed the correlation between other clinical factors and the immune profile, but found no statistically significant difference (Supplementary Figures 1B,C).
The association of stromal and immune scores with HCC prognosis was evaluated by dividing patients optimally into two groups based on their scores by using X-tile (see section “Materials and Methods” for details). We found that the high immune score and stromal score positively correlated with both OS and PFI (Figures 3A–D).
Figure 3. Kaplan–Meier (KM) survival curve of HCC patients based on their immune-stromal scores. Patients were classified into high immune-stromal scores groups and low immune/stromal scores groups by using X-tile. (A) The KM curve of overall survival (OS) time of high and low immune score group. (B) The KM curve of OS time of high and low stromal score group. (C) The KM curve of PFI time according to immune scores. (D) The KM curve of progression-free interval (PFI) time according to stromal scores. The survival outcomes of the two groups were compared by log-rank tests. P < 0.05 was statistically significant.
Comparison of Gene Expression Profile With Immune Scores and Stromal Scores in HCC
To identify the immune-related and stromal-related genes, differential analysis by using NetworkAnalyst was performed (Supplementary Table 3). The expression profiles of stromal and immune score–related DEGs are visualized, respectively, on the heatmaps (Figures 4A,B).
Figure 4. Expression profiles and biological functions of stromal and immune score–related DEGs. (A,B) Heatmaps showing expression profiles for selected stromal score (right) and immune score (left)–related DEGs (Log2 fold change ≥ 3 and adjusted P < 0.05) with unsupervised hierarchical clustering analyses, using the complete linkage method to measure distances between clusters. (C) Shows the commonly upregulated DEGs, and (D) shows the commonly downregulated DEGs. (E) The top 10 of biological processes GO terms (top), cellular component GO terms (middle), and molecular function GO terms (bottom); (F) KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis of microenvironment-related DEGs.
There were 797 shared DEGs overexpressed in both the stromal score and immune score groups (Figure 4C), and a total of 28 common DEGs were found to be underexpressed in both the stromal score and immune score groups (Figure 4D). Eight hundred twenty-five intersection genes were selected for further analysis (overlap zone in Figures 4C,D).
Using the “Goseq” and “clusterProfiler” R packages, 1,371 GO terms and 73 KEGG terms were indicated (Supplementary Table 4).
The results showed the top 10 biological processes GO terms, cellular component GO terms, and molecular function GO terms (Figure 4E). The correlation between the intersection genes and the top five biological processes is shown in Supplementary Figure 2A. The top 20 KEGG analysis showed that the intersection genes were associated with immune responses (Figure 4F).
Protein–Protein Interactions Among Intersection Genes
To better understand the interplay among the identified DEGs, we obtained PPI networks using the STRING tool. Using the MCODE software, we found modules in the network; the network was made up of eight modules, which included 408 nodes and 2,702 edges. We selected the top two significant modules for further analysis (Figure 5A and Supplementary Figure 2B).
Figure 5. Protein–protein interaction (PPI) network of microenvironment-related genes. (A) Module 1 is the top module in the PPI network. (B) GO analyses of module 1 (top 10 of biological processes GO terms). The color and thickness of edges reflect the combined score.
GO analyses of module 1 (Figure 5A) by ClueGo are shown in Figure 5B. Likewise, GO analyses of module 2 (Supplementary Figure 2B) by ClueGo are shown in Supplementary Figure 2C. The results demonstrated that module 1 was mainly enriched in regulation of dendritic cell apoptotic process, regulation of dendritic cell dendrite assembly, and positive regulation of T cell migration. Module 2 was mainly enriched in the regulation of phospholipase C activity, cellular response to interferon-γ (IFN-γ) and IFN-γ–mediated signaling pathway. Obviously, the top two modules were enriched for functional terms related to immune response processes, especially T cell responses.
Identification of Prognostic DEGs in HCC
To enrich for genes with the greatest prognostic values, we performed LASSO algorithm, and seven genes were identified (Supplementary Figure 3A). We also analyzed the association between the seven genes and OS using the Kaplan–Meier survival analysis. We found that the high levels of GDF10 (P = 0.0484) and MMP9 (P = 0.0143) negatively correlated with OS (Figure 6A).
Figure 6. Selection of microenvironment-related prognostic genes and the analysis of immune cell infiltration and immunoinhibitor. (A) Kaplan–Meier (KM) survival curve of GDF10 and MMP9. Patients were divided into two groups based on the median of gene expression. The survival outcomes of the two groups were compared by log-rank tests. P < 0.05 was statistically significant. (B) Correlation of microenvironment-related prognostic genes’ expression with immune infiltration level. (C) Relations between three kinds of immunoinhibitors and expression of MMP9. P < 0.05 was statistically significant, and partial correlation ≥0.3 indicates strong correlation. (D) The CARE score of MMP9 on CCLE, CGP, CTRP dataset. A positive score indicates a higher expression value to be associated with drug response.
Immune Cell Infiltration Analysis
To determine whether there is a correlation between tumor infiltration with immune cells and immune-related gene expression, the tumor infiltration with multiple immune cells was analyzed by TIMER 2.0 and other methods (Supplementary Table 5). Figure 6B shows the strong correlation between six types of immune cell infiltration and the expression of MMP9. The expression of MMP9 positively correlated with the infiltrating levels of B cells (partial correlation = 0.529, P = 3.05e-26), CD8+ T cells (partial correlation = 0.421, P = 4.13e-16), CD4+ T cells (partial correlation = 0.356, P = 9.68e-12), macrophages (partial correlation = 0.473, P = 2.12e-20), neutrophils (partial correlation = 0.34, P = 8.96e-11), and dendritic cells (partial correlation = 0.584, P = 1.72e-32). GDF10 expression was weakly associated with different immune cell infiltrates (Figure 6B).
We analyzed the correlation between MMP9 and immune checkpoints in liver cancer. MMP9 was found to be correlated with the expression of a series of immune checkpoints. Particularly, MMP9 was significantly correlated with PDCD1 (ρ = 0.576), PDCD1LG2 (ρ = 0.372), and CTLA4 (ρ = 0.672) (Figure 6C).
Besides, identifying reliable drug response biomarkers is a significant challenge in cancer research. We present CARE, a computational method that enables large-scale inference of response biomarkers and drug combinations for targeted therapies using compound screen data. High expression of MMP9 has been associated with better response to immunotherapies on CTRP dataset (Figure 6D).
Discussion
Prognosis prediction for liver cancer patients remains challenging for clinicians and investigators. Through a specific view of the microenvironment, this study provides a stromal-immune score–based gene signature to help answer this important clinical question.
Using the ESTIMATE algorithm, we revealed the correlation between the immune-stromal scores and the clinical HCC characteristics obtained from TCGA-CDR. The stromal and immune scores for tumor tissue were found to be positively associated with the clinicopathologic characteristics of the tumor and the patient’s prognosis. By analyzing the correlation between the immune scores and tumor recurrence, our data show that high-immune-score patients have a longer PFI and OS rates, indicating that the TME composition affects the clinical outcomes of HCC patients, which is consistent with previous studies (Haider et al., 2020).
Next, we analyzed 825 DEGs yielded from a comparison of high- versus low-immune-score (or stromal scores) groups and found that many of them were involved in the TME, specifically regulate T cell functions (Figure 4E). This is consistent with previous reports that the functions of immune cells and extracellular matrix molecules are interrelated in building TME in HCC (Lu et al., 2019; Yin et al., 2019). Moreover, we were able to construct two PPI modules (Figure 5 and Supplementary Figures 2B,C), the major of which were related to IFN-γ. We infer that these TME-associated genes might affect the development of HCC by affecting the T cell functions.
Finally, by using the LASSO algorithm (Supplementary Figure 3A), we identified seven TME-related genes. Of the seven genes identified, high levels of GDF10 and MMP9 showed a negative correlation to OS, which has been reported to be involved in carcinogenesis and the development of various cancers (Chang et al., 2017; Reggiani et al., 2017; Tekin et al., 2020a). We further correlated the degree of infiltration of six immune cell types with the expression of GDF10 and MMP9 by using TIMER algorithm. The expression of MMP9 was positively associated with the abundance of six immune in tumor tissues. It is worth reminding that our results did not contradict previous findings that high infiltration of CD8+ T cells indicated beneficial prognosis, but extended and enriched this conclusion. In the recent literature, tumor with higher CD8+ T cell infiltration, but T cell dysfunction and increased immune escape result in a poor prognosis (Hossain et al., 2020; Saka et al., 2020).
Prior studies have largely focused on MMPs’ ability to promote the invasion and metastasis of cancer cells (Nart et al., 2010; Chen et al., 2012), while evidence is mounting that MMPs are highly associated with the microenvironment of tumors and immune cells (Kessenbrock et al., 2010; Li et al., 2016). For example, MMP9-cleaved osteopontin fragments contribute to tumor immune escape by inducing the expansion of myeloid-derived suppressor cells (Shao et al., 2017). Macrophages secrete MMP9 to induce mesenchymal transition, which supports the tumor-promoting role of macrophage influx (Tekin et al., 2020b). Besides, MMP9 is associated with neutrophil migration (Koymans et al., 2016). Our study confirms the above conclusions and has found that MMP9 might associate with T cell dysfunction, despite high CD8+ cytotoxic T lymphocyte infiltration.
In addition, we also observed that high expression of MMP9 indicated higher levels of immune inhibitors (immune checkpoints), better response to immunotherapies, and poor survival in partial HCC patients, which was in line with our above analysis that some HCC patients with high CD8+ T cell infiltration but with dysfunction were immunosuppressed. And previously, inhibition of MMP9 could modulate immunosuppression in tumor (Melani et al., 2007). We also compared the prediction effect between the other factors, such as AFP (Supplementary Figure 1A) and programmed cell death protein 1 (PDCD1) (Supplementary Figure 3B), whereas AFP is not a good predictor of the abundance of immune invasion in HCC tissues, and PDCD1 is weakly correlated with the prognosis of HCC. Hence, MMP9 may be an effective biomarker to evaluate the immune status of patients and predict the effectiveness of immunotherapy before treatment. However, this conclusion will need to be confirmed by clinical trials in the future.
In summary, from comprehensively analyzing the correlation between microenvironmental and genetic factors of TCGA database applied by ESTIMATE algorithm-based immune and stromal scores, we identified MMP9 as a potential TME-related biomarker of prognostic and immunotherapy response. However, because of the lack of large sequenced HCC cohort and prospective clinical trials that have received immunotherapy, the effect of MMP9 expression on the efficiency of immunotherapy in HCC patients remains concerned.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.
Author Contributions
SL, GY, and LL: conception and design of study. SL: acquisition of data and drafting the manuscript. SL, XZ, LZ, EH, and YS: analysis and interpretation of data. XZ, GY, and LL: revising the manuscript critically for important intellectual content. SL, GY, LL, XZ, LZ, EH, and YS: approval of the version of the manuscript to be published. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National Natural Science Foundation of China (Grant Nos. 81773008 and 81672756), Guangdong Province Universities and Colleges Pearl River Scholar Funded Scheme (2015), Guangdong Natural Science Foundation (Grant No. 2017A030311023), Guangzhou Science and Technology Program projects (Grant No. 201804010044), Local Innovative and Research Teams Project of Guangdong Pearl River Talents Program (Grant No. 2017BT01S131), and Postdoctoral Research Foundation of China (Grant No. 2019M660208).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.625236/full#supplementary-material
Footnotes
- ^ https://xenabrowser.net
- ^ https://www.networkanalyst.ca/NetworkAnalyst/home.xhtml
- ^ http://www.ehbio.com/test/venn/#/
- ^ https://hiplot.com.cn
- ^ http://cis.hku.hk/TISIDB/
- ^ http://care.dfci.harvard.edu/
References
Affo, S., Yu, L. X., and Schwabe, R. F. (2017). The role of cancer-associated fibroblasts and fibrosis in liver cancer. Annu. Rev. Pathol. 12, 153–186. doi: 10.1146/annurev-pathol-052016-100322
Bader, G. D., and Hogue, C. W. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics 4:2. doi: 10.1186/1471-2105-4-2
Bardou, P., Mariette, J., Escudie, F., Djemiel, C., and Klopp, C. (2014). jvenn: an interactive Venn diagram viewer. BMC Bioinformatics 15:293. doi: 10.1186/1471-2105-15-293
Barry, A. E., Baldeosingh, R., Lamm, R., Patel, K., Zhang, K., Dominguez, D. A., et al. (2020). Hepatic stellate cells and Hepatocarcinogenesis. Front. Cell. Dev. Biol. 8:709. doi: 10.3389/fcell.2020.00709
Bindea, G., Mlecnik, B., Hackl, H., Charoentong, P., Tosolini, M., Kirilovsky, A., et al. (2009). ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 25, 1091–1093. doi: 10.1093/bioinformatics/btp101
Boyero, L., Sanchez-Gastaldo, A., Alonso, M., Noguera-Ucles, J. F., Molina-Pinelo, S., and Bernabe-Caro, R. (2020). Primary and acquired resistance to immunotherapy in lung cancer: unveiling the mechanisms underlying of immune checkpoint blockade therapy. Cancers (Basel) 12:3729. doi: 10.3390/cancers12123729
Bruix, J., Reig, M., and Sherman, M. (2016). Evidence-based diagnosis, staging, and treatment of patients with hepatocellular carcinoma. Gastroenterology 150, 835–853. doi: 10.1053/j.gastro.2015.12.041
Camp, R. L., Dolled-Filhart, M., and Rimm, D. L. (2004). X-tile: a new bio-informatics tool for biomarker assessment and outcome-based cut-point optimization. Clin. Cancer Res. 10, 7252–7259. doi: 10.1158/1078-0432.CCR-04-0713
Chang, Y. C., Chan, Y. C., Chang, W. M., Lin, Y. F., Yang, C. J., Su, C. Y., et al. (2017). Feedback regulation of ALDOA activates the HIF-1alpha/MMP9 axis to promote lung cancer progression. Cancer Lett. 403, 28–36. doi: 10.1016/j.canlet.2017.06.001
Chen, R., Cui, J., Xu, C., Xue, T., Guo, K., Gao, D., et al. (2012). The significance of MMP-9 over MMP-2 in HCC invasiveness and recurrence of hepatocellular carcinoma after curative resection. Ann. Surg. Oncol. 19 Suppl 3, S375–S384. doi: 10.1245/s10434-011-1836-7
Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33, 1–22.
Galle, P. R., Foerster, F., Kudo, M., Chan, S. L., Llovet, J. M., Qin, S., et al. (2019). Biology and significance of alpha-fetoprotein in hepatocellular carcinoma. Liver Int. 39, 2214–2229. doi: 10.1111/liv.14223
Haider, T., Sandha, K. K., Soni, V., and Gupta, P. N. (2020). Recent advances in tumor microenvironment associated therapeutic strategies and evaluation models. Mater. Sci. Eng. C. Mater. Biol. Appl. 116:111229. doi: 10.1016/j.msec.2020.111229
Hossain, M. A., Liu, G., Dai, B., Si, Y., Yang, Q., Wazir, J., et al. (2020). Reinvigorating exhausted CD8(+) cytotoxic T lymphocytes in the tumor microenvironment and current strategies in cancer immunotherapy. Med. Res. Rev. 41, 156–201. doi: 10.1002/med.21727
Jiang, P., Lee, W., Li, X., Johnson, C., Liu, J. S., Brown, M., et al. (2018). Genome-scale signatures of gene interaction from compound screens predict clinical efficacy of targeted cancer therapies. Cell Syst. 6, 343–354 e5. doi: 10.1016/j.cels.2018.01.009
Jin, M. Z., and Jin, W. L. (2020). The updated landscape of tumor microenvironment and drug repurposing. Signal Transduct. Target Ther. 5:166. doi: 10.1038/s41392-020-00280-x
Kanwal, F., and Singal, A. G. (2019). Surveillance for hepatocellular carcinoma: current best practice and future direction. Gastroenterology 157, 54–64. doi: 10.1053/j.gastro.2019.02.049
Kassambara, A., Kosinski, M., Biecek, P., and Fabian, S. (2019). survminer: Drawing Survival Curves Using ‘ggplot2’. R package version 0.4.6. Available online at: https://CRAN.R-project.org/package=survminer (accessed July 25, 2020).
Kessenbrock, K., Plaks, V., and Werb, Z. (2010). Matrix metalloproteinases: regulators of the tumor microenvironment. Cell 141, 52–67. doi: 10.1016/j.cell.2010.03.015
Korenberg, M. J. (2006). Applications of nonlinear system identification in molecular biology. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2006, 256–259. doi: 10.1109/IEMBS.2006.259455
Koymans, K. J., Bisschop, A., Vughs, M. M., van Kessel, K. P., de Haas, C. J., and van Strijp, J. A. (2016). Staphylococcal superantigen-like protein 1 and 5 (SSL1 & SSL5) limit neutrophil chemotaxis and migration through MMP-inhibition. Int. J. Mol. Sci. 17:1072. doi: 10.3390/ijms17071072
Li, M., Xing, S., Zhang, H., Shang, S., Li, X., Ren, B., et al. (2016). A matrix metalloproteinase inhibitor enhances anti-cytotoxic T lymphocyte antigen-4 antibody immunotherapy in breast cancer by reprogramming the tumor microenvironment. Oncol. Rep. 35, 1329–1339. doi: 10.3892/or.2016.4547
Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 77, e108–e110. doi: 10.1158/0008-5472.CAN-17-0307
Li, T., Fu, J., Zeng, Z., Cohen, D., Li, J., Chen, Q., et al. (2020). TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 48, W509–W514. doi: 10.1093/nar/gkaa407
Lu, C., Rong, D., Zhang, B., Zheng, W., Wang, X., Chen, Z., et al. (2019). Current perspectives on the immunosuppressive tumor microenvironment in hepatocellular carcinoma: challenges and opportunities. Mol. Cancer 18:130. doi: 10.1186/s12943-019-1047-6
Melani, C., Sangaletti, S., Barazzetta, F. M., Werb, Z., and Colombo, M. P. (2007). Amino-biphosphonate-mediated MMP-9 inhibition breaks the tumor-bone marrow axis responsible for myeloid-derived suppressor cell expansion and macrophage infiltration in tumor stroma. Cancer Res. 67, 11438–11446. doi: 10.1158/0008-5472.CAN-07-1882
Nart, D., Yaman, B., Yilmaz, F., Zeytunlu, M., Karasu, Z., and Kilic, M. (2010). Expression of matrix metalloproteinase-9 in predicting prognosis of hepatocellular carcinoma after liver transplantation. Liver Transpl. 16, 621–630. doi: 10.1002/lt.22028
Paul Blanche, J.-F. D., and Jacqmin-Gadda, H. (2013). Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat. Med. 32, 5381–5397. doi: 10.1002/sim.5958
Reggiani, F., Labanca, V., Mancuso, P., Rabascio, C., Talarico, G., Orecchioni, S., et al. (2017). Adipose progenitor cell secretion of GM-CSF and MMP9 promotes a stromal and immunological microenvironment that supports breast cancer progression. Cancer Res. 77, 5169–5182. doi: 10.1158/0008-5472.CAN-17-0914
Ru, B., Wong, C. N., Tong, Y., Zhong, J. Y., Zhong, S. S. W., Wu, W. C., et al. (2019). TISIDB: an integrated repository portal for tumor-immune system interactions. Bioinformatics 35, 4200–4202. doi: 10.1093/bioinformatics/btz210
Saka, D., Gokalp, M., Piyade, B., Cevik, N. C., Arik Sever, E., Unutmaz, D., et al. (2020). Mechanisms of T-Cell exhaustion in pancreatic cancer. Cancers (Basel) 12:2274. doi: 10.3390/cancers12082274
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Shao, L., Zhang, B., Wang, L., Wu, L., Kan, Q., and Fan, K. (2017). MMP-9-cleaved osteopontin isoform mediates tumor immune escape by inducing expansion of myeloid-derived suppressor cells. Biochem. Biophys. Res. Commun. 493, 1478–1484. doi: 10.1016/j.bbrc.2017.10.009
Son, J., Cho, J. W., Park, H. J., Moon, J., Park, S., Lee, H., et al. (2020). Tumor-infiltrating regulatory T cell accumulation in the tumor microenvironment is mediated by IL33/ST2 signaling. Cancer Immunol. Res. 8, 1393–1406. doi: 10.1158/2326-6066.CIR-19-0828
Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607–D613. doi: 10.1093/nar/gky1131
Tekin, C., Aberson, H. L., Waasdorp, C., Hooijer, G. K. J., de Boer, O. J., Dijk, F., et al. (2020a). Macrophage-secreted MMP9 induces mesenchymal transition in pancreatic cancer cells via PAR1 activation. Cell. Oncol. (Dordr.) 43, 1161–1174. doi: 10.1007/s13402-020-00549-x
Tekin, C., Aberson, H. L., Waasdorp, C., Hooijer, G. K. J., de Boer, O. J., Dijk, F., et al. (2020b). Macrophage-secreted MMP9 induces mesenchymal transition in pancreatic cancer cells via PAR1 activation. Cell. Oncol. (Dordr) 43, 1161–1174. doi: 10.1007/s13402-020-00549-x
Therneau, T. M. (2020). A Package for Survival Analysis in R. R Package Version 3.1-11. Available online at: https://CRAN.R-project.org/package=survival (accessed September 9, 2020).
Wang, J. B., Li, P., Liu, X. L., Zheng, Q. L., Ma, Y. B., Zhao, Y. J., et al. (2020). An immune checkpoint score system for prognostic evaluation and adjuvant chemotherapy selection in gastric cancer. Nat. Commun. 11:6352. doi: 10.1038/s41467-020-20260-7
Xia, J., Fjell, C. D., Mayer, M. L., Pena, O. M., Wishart, D. S., and Hancock, R. E. (2013a). INMEX–a web-based tool for integrative meta-analysis of expression data. Nucleic Acids Res. 41, W63–W70. doi: 10.1093/nar/gkt338
Xia, J., Lyle, N. H., Mayer, M. L., Pena, O. M., and Hancock, R. E. (2013b). INVEX–a web-based tool for integrative visualization of expression data. Bioinformatics 29, 3232–3234. doi: 10.1093/bioinformatics/btt562
Xu, D., Wang, Y., Liu, X., Zhou, K., Wu, J., Chen, J., et al. (2020). Development and clinical validation of a novel 9-gene prognostic model based on multi-omics in pancreatic adenocarcinoma. Pharmacol. Res. 164:105370. doi: 10.1016/j.phrs.2020.105370
Yang, J. D., Hainaut, P., Gores, G. J., Amadou, A., Plymoth, A., and Roberts, L. R. (2019). A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat. Rev. Gastroenterol. Hepatol. 16, 589–604. doi: 10.1038/s41575-019-0186-y
Yang, X., Chen, L., Liang, Y., Si, R., Jiang, Z., Ma, B., et al. (2018). Knockdown of alpha-fetoprotein expression inhibits HepG2 cell growth and induces apoptosis. J. Cancer Res. Ther. 14(Supplement), S634–S643. doi: 10.4103/0973-1482.180681
Yang, Y., Yang, Y., Yang, J., Zhao, X., and Wei, X. (2020). Tumor microenvironment in ovarian cancer: function and therapeutic strategy. Front. Cell. Dev. Biol. 8:758. doi: 10.3389/fcell.2020.00758
Yin, Z., Dong, C., Jiang, K., Xu, Z., Li, R., Guo, K., et al. (2019). Heterogeneity of cancer-associated fibroblasts and roles in the progression, prognosis, and therapy of hepatocellular carcinoma. J. Hematol. Oncol. 12:101. doi: 10.1186/s13045-019-0782-x
Yoshihara, K., Shahmoradgoli, M., Martinez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612
Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11:R14. doi: 10.1186/gb-2010-11-2-r14
Yu, G. (2019). enrichplot: Visualization of Functional Enrichment Result. R package version 1.6.1. Available online at: https://github.com/GuangchuangYu/enrichplot (accessed December 30, 2020).
Yu, G. (2020a). aplot: Decorate a ‘ggplot’ with Associated Information. R package Version 0.0.6. Available online at: https://CRAN.R-project.org/package=aplot (accessed September 3, 2020).
Yu, G. (2020b). Using ggtree to visualize data on tree-like structures. Curr. Protoc. Bioinformatics 69:e96. doi: 10.1002/cpbi.96
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118
Zhang, Y., Yang, M., Ng, D. M., Haleem, M., Yi, T., Hu, S., et al. (2020). Multi-omics data analyses construct TME and identify the immune-related prognosis signatures in human LUAD. Mol. Ther. Nucleic Acids 21, 860–873. doi: 10.1016/j.omtn.2020.07.024
Keywords: liver cancer, ESTIMATE, bioinformatics analysis, biomarker, tumor-microenvironment
Citation: Liu S, Yu G, Liu L, Zou X, Zhou L, Hu E and Song Y (2021) Identification of Prognostic Stromal-Immune Score–Based Genes in Hepatocellular Carcinoma Microenvironment. Front. Genet. 12:625236. doi: 10.3389/fgene.2021.625236
Received: 02 November 2020; Accepted: 06 January 2021;
Published: 11 February 2021.
Edited by:
Chunjie Jiang, University of Pennsylvania, United StatesReviewed by:
Yang Xie, Weill Cornell Medicine, United StatesMengqi Zhang, University of Pennsylvania, United States
Kuixi Zhu, University of Arizona, United States
Copyright © 2021 Liu, Yu, Liu, Zou, Zhou, Hu and Song. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Guangchuang Yu, gcyu1@smu.edu.cn; Li Liu, liuli@i.smu.edu.cn