Int J Med Sci 2020; 17(1):89-96. doi:10.7150/ijms.39056
Identification of gene expression profiles and immune cell infiltration signatures between low and high tumor mutation burden groups in bladder cancer
1. Department of Urology, Qilu Hospital of Shandong University, Jinan, China
2. Division of Gastroenterology, Tongji Hospital, Huazhong University of Science and Technology, Wuhan, China
*Zonglong Wu and Muru Wang contributed equally to this article.
Wu Z, Wang M, Liu Q, Liu Y, Zhu K, Chen L, Guo H, Li Y, Shi B. Identification of gene expression profiles and immune cell infiltration signatures between low and high tumor mutation burden groups in bladder cancer. Int J Med Sci 2020; 17(1):89-96. doi:10.7150/ijms.39056. Available from http://www.medsci.org/v17p0089.htm
Bladder cancer is one of the most commonly diagnosed tumors and is results from the accumulation of somatic mutations in the DNA. Tumor mutation burden (TMB) has been associated with cancer immunotherapeutic response. In this study, we attempted to explore the correlation between TMB and cancer prognosis. Identify the different expressed genes and immune cell infiltration signatures between low and high TMB group. Mutation data, gene expression profiles and clinical data were downloaded from The Cancer Genome Atlas (TCGA) database. Patients were divided into high and low TMB groups, allowing differentially expressed genes (DEGs) to be identified. Functional enrichment and protein-protein interaction (PPI) network analysis were used to identify the functions of the DEGs. And immune cell infiltration signatures were evaluated by CIBERSORT algorithm. These results shown that high TMB was significantly associated with prognosis. We obtained a list of TMB related genes which may influence the infiltrations of immune cells. We also found a higher proportion of CD8 T cells, CD4 T cells and NK cells in the high TMB group. Our data suggest that higher TMB tends to promote the infiltrations of T cells and NK cells and patients with higher TMB may achieve a more favorable prognosis in bladder cancer.
Keywords: tumor mutation burden, bladder cancer, gene expression profile, immune cell, TCGA database
Bladder cancer (BLCA) is the most common carcinoma in the urinary system and remains a major global medical problem despite the numerous new treatment options available, BLCA has a high recurrence rate which become a major economic burden on the health care systems [1, 2]. The cause of BLCA is unclear. Smoking, environmental factors, exposure to toxic industrial chemicals and gases, bladder inflammation are thought to be associated with bladder cancer . Despite improvements in clinical outcomes in recent years, there are still many issues to be aware of, such as patients with invasive tumors or metastatic have short survival times. Recently, immunotherapy target the programmed cell death 1 programmed cell death 1 (PD-1), PD-1 ligand (PD-L1), cytotoxic T-lymphocyte antigen-4 (CTLA-4) offers great promise for various cancer therapies [4-6].
Many studies have explored the link between immunotherapy response and tumor mutation burden (TMB) [7-10]. Mutations in tumor cells can be transcribed and translated, and may form new antigens can be identified and targeted by the immune system [11, 12]. In fact, not all mutations will generate immunogenic, only a few mutations can be recognized by T cells [13-15]. The more tumor mutations, the more antigens it may form. Higher TMB tends to form more new antigens, making tumors more immunogenic and improving clinical response to immunotherapy . So TMB can be used to estimate the new antigen load of a tumor.
The Cancer Genome Atlas (TCGA) has been established to map the genome variations of human cancers using genomic analysis techniques, providing a wealth of mutation and expression profile data. In this study, we identified different expression genes (DEGs) between high and low TMB groups and evaluate the relationship of immune cell infiltration signatures and TMB by using the information of patients with bladder cancer from the TCGA database.
Materials and Methods
Database and genomic analysis
The mutation data, gene expression profiles and clinical data of patients with bladder cancer were obtained from the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/) and the maftools package was used to analyze and summarized the mutation data. TMB was calculated from the tumor specific mutation genes. Gene expression data analysis was performed using the R software package, limma. A fold change of > 1.5 and false discovery rate (FDR) of < 0.05 were used as cutoffs to identify DEGs. Volcano plots and heat maps were generated using the ggplot2 and pheatmap packages respectively.
DEGs enrichment analysis and protein-protein interaction (PPI) network construction
Gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed using the clusterProfiler package . The STRING database  and cytoscape software  were used to retrieve and reconstruct a PPI network. Important nodes and subnetworks were predicted and explored using cytohubba, a cytoscape plugin , and the top 10 hub genes were selected from the results of each method.
Gene set enrichment analysis (GSEA) analysis
GSEA was used to further understand CXCL10-related pathways. The expression level of CXCL10 was used as the phenotype label and bladder cancer patients in the TCGA cohort were divided into two groups based on the median expression value of CXCL10. The collection of annotated gene sets of c2.cp.kegg.v6.2.symbols.gmt was chosen as the reference gene sets and GSEA version 3.0 software was used to analyze our data. FDR < 0.01 was used as the cut-off criteria.
Evaluation of tumor-infiltrating immune cells
CIBERSORT algorithm was used to calculate the fractions of infiltrating immune cells. CIBERSORT is an analytical tool that estimates specific cell types in a mixed cell population using gene expression data. And the algorithm was run using the 1000 permutations and LM22 signature . CIBERSORT method was used to quantify the fractions of immune cells in the bladder cancer samples. At a threshold of P < 0.05, the results of the inferred fractions of immune cell populations produced by CIBERSORT were considered accurate .
We divided patients into the high TMB or the low TMB groups with the median of mutation frequency as the threshold value. Kaplan‐Meier survival curves were obtained and compared by log‐rank tests. The associations of clinicopathologic characteristics and corresponding TMB were analyzed by one-way ANOVA followed by Tukey's multiple-comparison post-hoc test and unpaired two-tailed t test. The difference of infiltrating immune cells between high TMB group and low TMB group was assessed using the unpaired t test.
Primary genetic alterations in bladder cancer patients
In this study, the clinical information and the results of whole-exome sequencing of patients with BLCA were downloaded from the TCGA database. By using maftools, mutation data were analyzed and summarized. The mutations were further classified according to the variant effect predictor, among these mutations, missense mutations are the most common (Figure 1A). And the most common mutations type is SNP (Figure 1B). C > T transversion is the most common type of SNV in bladder cancer (Figure 1C). And the top 10 mutated genes are TTN, TP53, MUC16, KMT2D, ARID1A, KDM6A, SYNE1, PIK3CA, RB1, HMCN1(Figure 1D, E).
Correlation of TMB with prognosis, clinicopathological characteristics and tumor grades of BLCA patients
TMB was calculated as the number of nonsynonymous protein coding variants divided by the total sequenced genome length. Next, we divided the patients with bladder cancer into high and low TMB groups based on the median TMB. The clinicopathologic characteristics of the patients are shown in Table 1. Kaplan-Meier survival analysis revealed that patients with high TMB had a higher survival rate than those with low TMB (Figure 2A). Next, we analyzed the relationship between TMB and clinical stage, and the results showed that TMB had no relationship with the clinical stage (Figure 2B). And higher TMB level correlated with advanced tumor grades (Figure 2C).
Comparison of the gene expression profiles of patients in different TMB groups
Patients were divided into low and high TMB groups and their gene expression profiles were analyzed to identify DEGs with FDR < 0.05 and fold change of > 1.5. A total of 266 DEGs (89 up-regulated and 177 down-regulated) were identified in high TMB group (Figure 3A) and visualized using a heatmap (Figure 3B). And the list of DEGs is shown in supplementary Table S1.
Primary genetic alterations in bladder cancer patients. (A, B) Variant classification and type of genetic alterations in bladder cancer. (C). The SNV class of bladder cancer. (D, E) Top 10 mutant genes and mutation profile of bladder cancer.(Click on the image to enlarge.)
Correlation of TMB with prognosis, clinicopathological characteristics and tumor grades of BLCA patients. (A) Patients with BLCA were divided into two groups based on their TMB. As shown in the Kaplan‐ Meier survival curve, patients with high-TMB had a higher overall survival than those with low-TMB (hazard ratio [HR] 1.562; 95 % CI 1.14-2.14; P= 0.005 by log-rank test). (B) The TMB showed no statistically significant differences at different pathological stages (by one-way ANOVA followed by Tukey's multiple-comparison post-hoc test). (C) Higher TMB level correlated with advanced tumor grades (*, P<0.05; by unpaired two-tailed t test)(Click on the image to enlarge.)
Functional enrichment and PPI network analysis of differentially expressed genes
GO enrichment analysis was used to determine the functions of the 266 DEGs (Figure 4A). In BP category, “cell chemotaxis”, “lymphocyte chemotaxis” and “regulation of ion transmembrane transport” were enriched, which means the DEGs affects the consists of immune cells in tumor microenvironment. The enriched CC terms included “extracellular matrix”, “collagen-containing extracellular matrix”, “apical dendrite”, and the enriched MF terms included “chemokine activity”, “chemokine receptor binding”, “G protein-coupled receptor binding”. We also performed KEGG pathway enrichment analysis to determine the pathways most enriched for DEGs, which included “TGF-beta signaling pathway”, “chemokine signaling pathway” and “pathways in cancer” (Figure 4B). Next, we explored the relationships between the DEGs, The Search Tool for the Retrieval of Interacting Genes (STRING) database and Cytoscape software were used to construct a PPI network for the DEGs (Figure 5A). The important nodes and subnetworks of the PPI were predicted and explored using CytoHubba; the 10 most significant node genes were CXCL10, CXCL11, GNG7, CXCR2, AGT, ADCY5, CCL5, ADRA2A, S1PR1, GALR2. Next, we analyzed the pathway of high-expression samples of CXCL10 by GSEA analysis. The results showed that high CXCL10 expression samples were mainly enriched in natural killer cell mediated cytotoxicity, antigen processing and presentation, chemokine signaling pathway and T cell receptor signaling pathway (Figure 5B).
Correlation of TMB with immune signatures in bladder cancers
Previous studies have shown that the higher mutation burden in tumors tends to form more new antigens, making tumors to have higher immunogenicity . The GO enrichment analysis shows the DEGs involve in cell chemotaxis, so we analyzed the correlation of TMB with immune signatures in bladder cancers. CIBERSORT algorithm was used to calculated the fractions of infiltrating immune cells. At a threshold of P < 0.05, the results of the inferred fractions of immune cell populations produced by CIBERSORT were considered accurate. There are 90 patients in low TMB group and 104 patients in high TMB group with the P value < 0.05 and their inferred fractions of immune cell were considered accurate. The results shown that tumors with high TMB were significantly associated with high fractions of CD8 T cells, CD4 memory T cells, follicular helper T cells and resting NK cells. In low TMB group there is a higher fraction of mast cell. (Figure 6).
Comparison of the gene expression profiles of patients in different TMB groups. (A, B) The volcano plot and heatmap show the 266 genes (89 up-regulated and 177 down-regulated) identified based on the TMB.(Click on the image to enlarge.)
Functional enrichment of differentially expressed genes. (A) Biological process, cellular component, and molecular function terms for the DEGs. (B) KEGG pathways enriched for the DEGs.(Click on the image to enlarge.)
PPI network analysis of differentially expressed genes. (A) Protein-protein interaction networks of the DEGs. (B) CXCL10 correlated enrichment gene analysis with GSEA.(Click on the image to enlarge.)
Cancer is a genetic disease and results from the accumulation of somatic mutations in the DNA . Genetic changes in tumors include nonsynonymous mutations, synonymous mutations, insertions or deletions, and copy number gains and losses. And nonsynonymous mutations mainly comprised of missense mutations (point mutations that change the amino acid codon). In different tumor types and individual tumors, there is obvious difference in the frequency of each type of these genetic alterations . The genetic changes increase the tumors' immunogenic. To avoid detection and killing by the host immune system, tumors often upregulate immune checkpoints. Recently, the overall survival (OS) rates of bladder cancer have increased with the therapy of immune checkpoint . TMB can be used to predict the efficacy of immune checkpoint blockade therapy, and can been seen as a useful biomarker to identify the patients who will benefit from immunotherapy [9, 10, 24].
High TMB cases can be seen in almost every type of cancer . In different tumor types, melanomas have the highest levels of TMB followed by non-small-cell lung carcinoma and other squamous carcinomas. Leukemias and pediatric tumors usually have the lowest levels of TMB . And TMB also has a significant difference in the same cancer type. A high TMB probably reflects the presence of mutation-associated neoantigens, with consequent increased lymphocyte infiltration in the tumor microenvironment. This phenomenon has been observed in other tumors [26, 27]. The tumor microenvironment consists of immune cells, mesenchymal cells, endothelial cells, extracellular matrix (ECM) molecules, and inflammatory mediators. BLCA is an immunosensitive tumor which is infiltrated by tumor-infiltrating immune cells (TIICs) including T cells, macrophages, dendritic cells, neutrophils and mast cells [28-30]. Studies have shown that the tumor microenvironment affects the gene expression of tumor tissues and the patient outcome, and therefore, has a diagnostic and prognostic value for BLCA . In present study, we found that TMB affects the prognosis of bladder cancer. High TMB may reflect the presence of new antigens, thereby increasing tumor-infiltrating immune cells in the tumor microenvironment which closely related to the effectiveness of targeted drugs and clinical outcomes. We also identified 266 DEGs (89 up-regulated and 177 down-regulated) between low and high TMB groups. Among the DEGs, many genes involve immune response and chemokine signaling pathway. GO term and KEGG analysis revealed that the DEGs affects the cell chemotaxis, intercellular signaling, ion transport, and the formation of extracellular matrix. These data indicated that TMB is closely related with the tumor microenvironment and these TMB-related genes cause the changes of tumor microenvironment. CXCL10 is the most significant node genes Among all DEGs and high CXCL10 expression samples were mainly enriched in natural killer cell mediated cytotoxicity, antigen processing and presentation, chemokine signaling pathway and T cell receptor signaling pathway. Chemokines and cytokines are well known to guide macrophages, T-cells and other immune cells to the tumor microenvironment and influence the outcome of the patients . CXCL10 may play an important role in regulating immune cell migration, differentiation, and activation in bladder cancer.
Correlation of TMB with immune signatures in bladder cancers. 22 types of adaptive and innate immune cells in high and low TMB groups. (*, P<0.05; **, P<0.01; ***, P<0.001; by unpaired two-tailed t test).(Click on the image to enlarge.)
Immune-cell infiltration is a characteristic of cancer, and many cancers have a complex chemokine network that affect the extent and phenotype of this infiltrate, as well as tumor cell growth, survival and migration . In this study, we found high TMB group has higher fractions of CD8 T cells, CD4 T cells and NK cell and in low TMB group; there is a higher fraction of mast cell. These data indicated that TMB can affect the immune cell infiltration signatures and high TMB attracted effector cells of the immune system. Tissue resident memory T cells are a key factor in making tumors dormant; hence, it is essential to establish a cancer-immune system balance . Under hypoxic conditions, CD8 T cells can differentiate into lytic effector cells, increase the expression of interferon gamma (IFNγ), Fas ligand (FASL), granule B (GZMB), and inhibit tumor cell proliferation [35, 36]. High infiltration of T lymphocytes in tumors is positively correlated with the survival rate of patients with bladder cancer . NK cell typically account for 5-15% of peripheral blood lymphocytes and respond to their targets without prior antigen sensitization . NK cells can recognize bladder tumor cells and their activity is important to against bladder tumor cells . Mast cells may contribute to tumor angiogenesis and play an important role in the growth of tumors . These data indicated that TMB is closely related with the immune microenvironment. Mast cells in the low TMB group may promote tumor growth and metastasis. High TMB tends to cause the chemotaxis of immune cells in BLCA and the crosstalk between these cells play an important role in the growth of tumors.
In summary, our data implicate that higher-TMB patients could gain a more favorable prognosis in bladder cancer. We also obtained a list of TMB related genes which may influence the infiltrations of immune cells. Our data provide insights into the correlation between TMB and immune cell infiltration signatures in bladder cancer and may be helpful for the exploration of the role of TMB in BLCA.
Supplementary figures and tables.
This study was supported by the Funds of the National Natural Science Foundation of China (Grant No. 81470987 and NO. 81900637) and the Tai Shan Scholar Foundation to Benkang Shi.
The authors have declared that no competing interest exists.
1. Antoni S, Ferlay J, Soerjomataram I, Znaor A, Jemal A, Bray F. Bladder Cancer Incidence and Mortality: A Global Overview and Recent Trends. Eur Urol. 2017;71(1):96-108
2. Dobruch J, Daneshmand S, Fisch M. et al. Gender and Bladder Cancer: A Collaborative Review of Etiology, Biology, and Outcomes. Eur Urol. 2016;69(2):300-10
3. Janković S, Radosavljević V. Risk factors for bladder cancer. Tumori. 2007;93(1):4-12
4. Zhang J, Bu X, Wang H. et al. Cyclin D-CDK4 kinase destabilizes PD-L1 via cullin 3-SPOP to control cancer immune surveillance. Nature. 2018;553(7686):91-95
5. Kim HS, Seo HK. Immune checkpoint inhibitors for urothelial carcinoma. Investig Clin Urol. 2018;59(5):285-296
6. Powles T, Morrison L. Biomarker challenges for immune checkpoint inhibitors in urothelial carcinoma. Nat Rev Urol. 2018;15(10):585-587
7. Van Allen EM, Miao D, Schilling B. et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science. 2015;350(6257):207-211
8. June CH, O'Connor RS, Kawalekar OU, Ghassemi S, Milone MC. CAR T cell immunotherapy for human cancer. Science. 2018;359(6382):1361-1365
9. Rizvi NA, Hellmann MD, Snyder A. et al. Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science. 2015;348(6230):124-8
10. Hugo W, Zaretsky JM, Sun L. et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell. 2017;168(3):542
11. Matsushita H, Vesely MD, Koboldt DC. et al. Cancer exome analysis reveals a T-cell-dependent mechanism of cancer immunoediting. Nature. 2012;482(7385):400-4
12. Riaz N, Morris L, Havel JJ, Makarov V, Desrichard A, Chan TA. The role of neoantigens in response to immune checkpoint blockade. Int Immunol. 2016;28(8):411-9
13. Coulie PG, Van den Eynde BJ, van der Bruggen P, Boon T. Tumour antigens recognized by T lymphocytes: at the core of cancer immunotherapy. Nat Rev Cancer. 2014;14(2):135-46
14. Carreno BM, Magrini V, Becker-Hapak M. et al. Cancer immunotherapy. A dendritic cell vaccine increases the breadth and diversity of melanoma neoantigen-specific T cells. Science. 2015;348(6236):803-8
15. Snyder A, Chan TA. Immunogenic peptide discovery in cancer genomes. Curr Opin Genet Dev. 2015;30:7-16
16. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284-287
17. Szklarczyk D, Franceschini A, Wyder S. et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43(Database issue):D447-52
18. Shannon P, Markiel A, Ozier O. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498-504
19. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4):S11
20. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol. 2018;1711:243-259
21. Chan TA, Yarchoan M, Jaffee E. et al. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol. 2019;30(1):44-56
22. Lawrence MS, Stojanov P, Polak P. et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013;499(7457):214-218
23. Chang SS. Re: Atezolizumab in Patients with Locally Advanced and Metastatic Urothelial Carcinoma Who Have Progressed following Treatment with Platinum-Based Chemotherapy: A Single-Arm, Multicentre, Phase 2 Trial. J Urol. 2016;196(6):1637-1638
24. Chan TA, Wolchok JD, Snyder A. Genetic Basis for Clinical Response to CTLA-4 Blockade in Melanoma. N Engl J Med. 2015;373(20):1984
25. Chalmers ZR, Connelly CF, Fabrizio D. et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 2017;9(1):34
26. Innocenti F, Ou FS, Qu X. et al. Mutational Analysis of Patients With Colorectal Cancer in CALGB/SWOG 80405 Identifies New Roles of Microsatellite Instability and Tumor Mutational Burden for Patient Outcome. J Clin Oncol. 2019;37(14):1217-1227
27. Wang X, Li M. Correlate tumor mutation burden with immune signatures in human cancers. BMC Immunol. 2019;20(1):4
28. Schwamborn K, Ammann JU, Knüchel R. et al. Multicentric analytical comparability study of programmed death-ligand 1 expression on tumor-infiltrating immune cells and tumor cells in urothelial bladder cancer using four clinically developed immunohistochemistry assays. Virchows Arch. 2019
29. Stone L. Bladder cancer: Mastering the immune microenvironment. Nat Rev Urol. 2017;14(11):639
30. Sweis RF, Spranger S, Bao R. et al. Molecular Drivers of the Non-T-cell-Inflamed Tumor Microenvironment in Urothelial Bladder Cancer. Cancer Immunol Res. 2016;4(7):563-568
31. Efstathiou JA, Mouw KW, Gibb EA. et al. Impact of Immune and Stromal Infiltration on Outcomes Following Bladder-Sparing Trimodality Therapy for Muscle-Invasive Bladder Cancer. Eur Urol. 2019;76(1):59-68
32. Galon J, Angell HK, Bedognetti D, Marincola FM. The continuum of cancer immunosurveillance: prognostic, predictive, and mechanistic signatures. Immunity. 2013;39(1):11-26
33. Balkwill F. Chemokine biology in cancer. Semin Immunol. 2003;15(1):49-55
34. Gabriel SS, Kallies A. Tissue-resident memory T cells keep cancer dormant. Cell Res. 2019;29(5):341-342
35. de Silly R V, Ducimetière L, Yacoub MC, Dietrich PY, Derouazi M, Walker PR. Phenotypic switch of CD8(+) T cells reactivated under hypoxia toward IL-10 secreting, poorly proliferative effector cells. Eur J Immunol. 2015;45(8):2263-2275
36. Xu Y, Chaudhury A, Zhang M. et al. Glycolysis determines dichotomous regulation of T cell subsets in hypoxia. J Clin Invest. 2016;126(7):2678-2688
37. Winerdal ME, Marits P, Winerdal M. et al. FOXP3 and survival in urinary bladder cancer. BJU Int. 2011;108(10):1672-8
38. Vivier E, Raulet DH, Moretta A. et al. Innate or adaptive immunity? The example of natural killer cells. Science. 2011;331(6013):44-9
39. García-Cuesta EM, López-Cobo S, Álvarez-Maestro M. et al. NKG2D is a Key Receptor for Recognition of Bladder Cancer Cells by IL-2-Activated NK Cells and BCG Promotes NK Cell Activation. Front Immunol. 2015;6:284
40. Sari A, Calli A, Cakalagaoglu F, Altınboga AA, Bal K. Association of mast cells with microvessel density in urothelial carcinomas of the urinary bladder. Ann Diagn Pathol. 2012;16(1):1-6
Corresponding authors: Yan Li, Department of Urology, Qilu Hospital of Shandong University, 107 Wenhuaxi Road, Jinan, 250012, China, E-mail: yanliedu.cn. Benkang Shi, Department of Urology, Qilu Hospital of Shandong University, 107 Wenhuaxi Road, Jinan, 250012, China, E-mail: bkang68edu.cn.