Establishment of a prognostic model related to tregs and natural killer cells infiltration in bladder cancer

Abstract

BACKGROUND

Regulatory T cells (Tregs) and natural killer (NK) cells play an essential role in the development of bladder urothelial carcinoma (BUC).

AIM

To construct a prognosis-related model to judge the prognosis of patients with bladder cancer, meanwhile, predict the sensitivity of patients to chemotherapy and immunotherapy.

METHODS

Bladder cancer information data was obtained from The Cancer Genome Atlas and GSE32894. The CIBERSORT was used to calculate the immune score of each sample. Weighted gene co-expression network analysis was used to find genes that will have the same or similar expression patterns. Subsequently, multivariate cox regression and lasso regression was used to further screen prognosis-related genes. The prrophetic package was used to predict phenotype from gene expression data, drug sensitivity of external cell line and predict clinical data.

RESULTS

The stage and risk scores are independent prognostic factors in patients with BUC. Mutations in FGFR3 lead to an increase in Tregs percolation and affect the prognosis of the tumor, and additionally, EMP1, TCHH and CNTNAP3B in the model are mainly positively correlated with the expression of immune checkpoints, while CMTM8, SORT1 and IQSEC1 are negatively correlated with immune checkpoints and the high-risk group had higher sensitivity to chemotherapy drugs.

CONCLUSION

Prognosis-related models of bladder tumor patients, based on Treg and NK cell percolation in tumor tissue. In addition to judging the prognosis of patients with bladder cancer, it can also predict the sensitivity of patients to chemotherapy and immunotherapy. At the same time, patients were divided into high and low risk groups based on this model, and differences in genetic mutations were found between the high and low risk groups.

Key Words: Natural killer cells; Tregs; Bladder cancer; Weighted gene coexpression network analysis; Bladder cancer treatment; Immunotherapy; Computational molecular biology

Core Tip: Tregs are thought to be connected to tumor cells evading the immune system, which is linked to cancer patients' poor prognoses. Natural killer (NK) cells control different immune responses and show antitumor cytotoxicity without being previously sensitized. We determined the immunological scores of several cell types in bladder urothelial carcinoma, discovered a gene set that was favorably connected with the Tregs score and negatively correlated with the NK cells score, and built a model that was related to prognosis. The model can predict the sensitivity of patients to chemotherapy and immunotherapy in addition to their prognosis for bladder cancers.

INTRODUCTION

Bladder urothelial carcinoma (BUC) is a typical malignant growth of the urinary tract that has a high likelihood of return and death[1]. Immunotherapy has lately attracted more research attention than traditional surgery methods and chemotherapy[2]. Understanding the changes in the tumor microenvironment and finding potential immunotherapy strategies have received a lot of attention from scientific specialists. Thus, it is crucial to research how different immune cells infiltrate the tumor's inflammatory environment.

The start and progression of BUC are influenced by three processes: Immune monitoring, immunological balance, and immune escape. T cells (Tregs) are currently believed to be implicated in tumor cell immune evasion; Tregs are frequently found in tumors and are typically identified by the expression of CD4⁺ and Foxp3⁺, which is linked to the bad outcome of cancer patients[3]. The grade and stage of BUC were considerably correlated with the proportion of Foxp3⁺ Tregs[4]. Inhibiting Tregs can improve the treatment of BUC patients. Natural killer (NK) cells are mainly related to killing infected microorganisms and malignant transformed allogeneic and autologous cells. NK cells exhibit antitumor cytotoxicity without prior sensitization and production of cytokines and chemokines that regulate various immune responses[5]. NK cells patrol the body and are recruited to tumor sites to destroy malignant cells[6].

Based on the infiltration of Treg cells and NK cells into the tumor tissue, we built a prognosis-related model of bladder tumor patients for this investigation. It can predict the susceptibility of patients to chemotherapy and immunotherapy in addition to assessing the prognosis of bladder cancer patients.

MATERIALS AND METHODS

Downloading and processing of data

The validated database is derived from GSE32894. The platform of the data set is GPL6947. In the Gene Expression Omnibus (GEO) database, the end-point event is determined according to the disease-specific survival. GSE32894 included 308 patients with BUC. The data source of the training data set is the transcriptome data of The Cancer Genome Atlas (TCGA's) bladder cancer, including 19 patients' normal samples and 411 patients' cancer samples. The Data Category and Workflow Type of TCGA data are transcriptome profiling and Fragments per Kilobase Million (FPKM). Then, the TranscriptsPerKilobase of exonmodel per Million mapped reads (TPM) method was used to standardize the TCGA data for model verification. The mutation data were obtained from the simple nucleoside variation data of 409 patients with bladder urothelial carcinoma in TCGA database.

Calculation of immune cell score

Each sample's immunological score was determined using "CIBERSORT". 22 immune cells in total were scored, and "perm = 1000, qn = true" was specified as the option in "CIBERSORT".

Weighted gene coexpression network analysis

Finding genes that will have the same or comparable expression patterns may be done using the useful method known as weighted gene coexpression network analysis (WGCNA). These genes could have comparable physiological effects. WGCNA differs from other straightforward clustering methods in this way (such as based on Euclidean distance). It is a clustering technique with biologically relevant applications. In this work, we employed WGCNA to categorize the bladder cancer expression profile data into different groups.

Screening of prognostic related genes

Then, the prognosis-related diagnostic genes in the TCGA database were screened using the univariate cox regression model. Multivariate cox regression and lasso regression were utilized to further screen prognosis-related genes based on the findings of the univariate cox screening of the TCGA database and GEO data. We determined the riskcore for each sample using the multivariate Cox model (riskcore = gene expression * coef).

Drug sensitivity analysis

Prrophetic package is an ancient R package. Its main purpose is to predict phenotype from gene expression data (predict clinical results using CGP cell line data of cancer genome project), predict drug sensitivity of external cell line, and predict clinical data.

Statistical methods

The clinical data were compared using the chi square test. The correlation between different values was determined using Pearson correlation analysis. The following is the study's flowchart (Figure 1).

Figure 1 The flowchart of our study. TCGA: The Cancer Genome Atlas; NK: Natural killer; WGCNA: Weighted gene coexpression network analysis.

Waterfall diagram of gene mutation

Maftools package is used to analyze the overall situation of gene mutation between high-risk and low-risk groups. The top 20 genes with the highest frequency of mutations in both high-risk and low-risk groups are displayed in a waterfall chart.

The Cancer Immunome Atlas

The Cancer Immunome Atlas (TCIA) is developed based on TCGA data. The difference is that TCIA only provides immune data analysis of 20 cancer species. The immunophenoscore (IPS) of each sample in the high and low risk categories, as well as the difference between them, were determined using TCIA. The reactivity of CTLA-4 and PD-1 can be accurately predicted by IPS.

RESULTS

Immune cell score and WGCNA

In the TCGA database, we first assigned a number to each sample's immune cells (Figure 2A). The total genes were divided into 9 gene sets based on the immune cell scores of each sample using the WGCNA technique, and the correlation between these 9 gene sets and the immune cell scores of each type was determined (Figure 2B-E). Then we chose a group of genes (module yellow) that had a positive correlation with the score for Regulatory T cells (Treg) (r = 0.45) and a negative correlation with the score for NK cell activation (r = -0.33). According to earlier studies, Tregs inhibit long-lasting immune reactions to viruses, tumors, and self-antigens; the infiltrate tumor tissue and are linked to a bad prognosis in cancer patients; the quantity of Treg was negatively correlated with recurrence-free survival in BUC and substantially increased in peripheral blood and tumor tissue[7]. The grade and stage of BUC were significantly associated with the proportion of Foxp3+ Tregs[8]. The most significant component of natural immunity is played by NK cells, which can directly destroy target cells like cancerous or virus-infected ones. They serve as the first line of protection against tumors and the first obstacle of human antibodies. From this perspective, the gene collection we chose is also consistent with the findings of earlier studies. More significantly, prior research has shown that tumor-infiltrating Treg cells can cause NK cells, B cells, Dendritic cells (DC) and cytotoxic T cells in the Tumor microenvironment (TME) to undergo death[9]. This strengthens the logic of our study concepts even more.

Figure 2 Immune infiltration status of TCGA database samples and gene sets related to immune cell infiltration. A: Heatmap of the distribution of immune cell infiltration in each sample; B: The scale independence and the mean connectivity; C: Hierarchical clustering was used to obtain coexpressive gene modules; Different gene modules are assigned to different colors; D: Clustering diagram of each gene module; E: Clustering diagram of each gene module.

Screening of prognosis related genes and construction of model

By using multivariate Cox regression analysis, the prognostic-related genes in the yellow gene group were chosen (Figure 3A). The alleles with high association were then eliminated using lasso regression (Figure 3B). Lastly, the prognostic association model was built using multivariate Cox regression. 13 genes in total (CHAT, EMP1, ADSl, TRIM16 L, CMTM8, DNASE2B, NTN4, SORT1, IQSEC2, TCHH, CNTNAP3B, ESD, NPSR1) were used to build this model (Figure 4A and B). The riskscore of each sample was then determined using the multivariate Cox model (riskscore = gene expression * coef). Samples that scored above the median for hazards were placed in the high-risk category, while those that scored below the median were placed in the low-risk group. Then it was looked at how long the high-risk group survived compared to the low-risk group. Clearly, the prognosis between the high-risk group and the low-risk group is extremely different (Figure 3C). In comparison to the high-risk group, the prognosis for the low-risk group was noticeably better. The GEO database has verified the model's efficacy in determining prognosis (Figure 4D). In the supporting materials, a survival analysis of individual genes will be presented. Furthermore, we determined the association (r > 0.4 and P value > 0.05 were deemed significant) between each gene's expression and immune cell infiltration in the model. The prognosis of patients was then predicted using Nomo map (Figure 2E; stage: "1","2","3","4" represent Stage I, Stage II, Stage III, and Stage IV, respectively; grade: "1" and "2"represent G1 and G2, respectively; gender: "0" and "1" represent females and males, respectively).

Figure 3 Screening of survival related gene sets (based on immune infiltration). A: Screening of prognostic genes by univariate Cox; B: Least absolute shrinkage and selection operator (LASSO) regression model; C: Kaplan Meier curve shows that in The Cancer Immunome Atlas (TCGA) database, the prognosis of patients between high and low risk groups is different (Training group); D: Kaplan Meier curve shows the difference in prognosis between high and low risk groups in geo database (Test group); E: Nomograms were used to measure patient survival in The Cancer Immunome Atlas database.

Figure 4 Differences in prognosis between high-risk and low-risk patients. A: Heat map of gene expression differences between high and low risk groups in The Cancer Genome Atlas database; B: Heat map of gene expression differences between high and low risk groups in Gene Expression Omnibus database; E: The ROC curve of the risk score and clinicopathological variables; F: Distribution heat map of model genes and clinicopathological variables; G: The ROC curve of the risk score and clinicopathological variables; H: Distribution heat map of model genes and clinicopathological variables.

Model independence and diagnostic efficiency

Cox regression analysis was conducted to determine whether our model could predict the prognosis of BUC patients independently of other clinical data. Stage and risk score were shown to be strongly correlated with the prognosis of BUC patients by univariate Cox regression analysis. Stage and risk score were identified as independent predictive variables in multivariate Cox regression analysis of BUC patients (Figure 4E and F). Furthermore, the scatter plot demonstrates that patients who scored highly for hazards have much poorer prognoses (Figure 4C and D). The receiver operating characteristic curve (AUC at 1 year: 0.699, AUC at 3 years: 0.745, AUC at 5 years: 0.753) further proves that the model can more accurately predict the prognosis of patients (Figure4G). We employed chi-square test to evaluate the difference in clinical data between high and low risk groups ("***" = P0.001, "**" = P0.01, "*" = P0.05) in order to further investigate the difference between high and low risk groups. Patients in the high-risk group underwent much more distant metastases and had significantly later clinical stages. (Figure 4H).

Differences in genetic mutations between high and low risk

The variations in genetic mutations between high and low risk were then looked at. According to the findings, the prevalence of FGFR3 mutations varied considerably between the high-risk and low-risk groups (Figure 5A and B). One of the frequent mutations in bladder cancer is the FGFR3 mutation, which is also one of the key elements in carcinogenesis. Moreover, FGFR3-targeting medications have been utilized to treat BUC patients. It will be fascinating to explore if the FGFR3 mutation would boost Treg infiltration and have an impact on the prognosis of tumors.

Figure 5 Different treatment-related characteristics of patients in the high-risk and low-risk groups. A and B: Significant differences in mutation frequency of gene between the high risk group and the low risk group; C: Correlation between genes, riskscore and immune checkpoints in the model; D: The Estimation of the immunophenoscore in Immunotherapy Response.

The difference of IPS (immunophenoscore) and gene mutation between high-risk and low-risk groups, as well as the difference between the model and each immune checkpoint

The model and each immune checkpoint differ in that EMP1, TCHH, and CNTNAP3B are mostly favorably connected with the expression of immunological checkpoints, whereas CMTM8, SORT1, and IQSEC1 are adversely linked with immune checkpoints (Figure 5C). Moreover, the IPS between the high-risk and low-risk groups varied (Figure 5D). It is evident that the immune checkpoint inhibition may work better for the low-risk population.

Prediction of drug sensitivity

To distinguish between high and low risk groups in terms of medication sensitivity, we employed the prrophetic package (R package). Interestingly, chemotherapeutic medication sensitivity was greater in the high-risk group. These chemotherapy medications include vinblastine, cisplatin, docetaxel, mitomycin, and adriamycin (Figure 6). These are all chemotherapeutic medications that are frequently used to treat bladder cancer.

Figure 6 Differences in IC50 values of different drugs between high and low risk groups. A: IC50 values of cisplatin between high and low risk groups; B: IC50 values of Docetaxel between high and low risk groups; C: IC50 values of vinblastine between high and low risk groups; D: IC50 values of Doxorubicin between high and low risk groups; E: IC50 values of Mitormycin. C between high and low risk groups.

DISCUSSION

Tregs, an immunosuppressive subset of CD4⁺ T cells, play a crucial role in maintaining self-tolerance, avoiding autoimmune disorders and transplant rejection, it has been proposed as candidates of bio-markers because of their capability to control alloimmune responses[10]. In tumor immunity, Tregs compromise immune surveillance against cancer in healthy individuals and damage the antitumor immune response in tumor-bearing hosts[11]. Among different types of tumor infiltrating lymphocytes, the infiltration of cytotoxic CD8⁺T cells plays a positive role in the prognosis of tumor patients and the clinical response of immunotherapy[12,13], but the infiltration of tumor infiltrating lymphocytes called Treg commonly indicates poor prognosis and poor effect of immunotherapy. Tregs have been reported to increase infiltration in various tumor tissues, which is related to tumor immune escape.

Tregs can secrete anti-inflammatory mediators including Interleukin (IL)-10,Transforming factor β (TGF-β), and IL-35, which can suppress the immune system. TGF-β lowers the lethality of NK cells and induces them to transform into type 1 innate lymphocytes, which has further effects on tumor development and metastasis[14].

Due to their quick recognition and effective killing of tumor cells, NK cells, which are a crucial component of innate immunity, serve as the first line of defense against tumors[15]. Meanwhile, NK cells play a central role in cancer immune surveillance through express up to different activating inhibitory killer-cell immunoglobulin-like receptors[16]. The NK cells float throughout the body. NK cells are recruited to the cancer site when the cell becomes cancerous. It plays a positive role in the prognosis of tumor patients, as has been demonstrated in several cancer syndromes[17]. In addition to killing tumor cells directly, NK cells can also promote the recruitment of cdc1s to TME by producing CCL5, Xcl1/2 and FLT3LG. Cdc1s is a DC subgroup, which is responsible for cross presentation of tumor antigens to CD8⁺T cells, indicating that NK cells play a key role in enhancing anti-tumor CD8⁺T cell response[18,19].

The number of Treg in bladder cancer tissue is significantly higher than that in non-cancerous tissue; Fattahi et al[8] proved that conclusion through immunohistochemical.

EMP1 has been reported to be a prooncogene in many cancers, such as bladder cancer, childhood leukemia, non-small cell lung cancer and glioma. Miao et al[20] reported that EMP1 promotes tumor progression by targeting C-MYC. In BLCA, there was a moderate positive correlation between the number of Tregs infiltrating and EMP1 expression and this result was calculated[21,22]. Taha-Mehlitz et al[23] found that ADSL affects mitochondrial function by altering the Tricarboxylic acid cycle and mitochondrial respiratory injury, thus playing its role in promoting the progression of colon cancer. Gao et al[24] showed that the expression of CMTM8 was down regulated in bladder cancer and existed as a tumor suppressor gene. Overexpression of CMTM8 could improve the sensitivity of bladder cancer cell lines to epirubicin. Moreover, the expression of CMTM8 has been reported to be decreased in a variety of tumors, such as liver cancer, colon, lung adenocarcinoma and so on. At the same time, CMTM8 is also associated with EMT (epithelial mesenchymal transformation) and immune infiltration[25,26]. The conclusions drawn from this article contradict our results. DNASE2B may be associated with chemotherapy sensitivity of gastric cancer[27]. Nerve guidance factor 4 (Netrin-4, NTN4) is a ligand of neo1 (a transmembrane receptor). Its expression in gastric cancer, breast cancer and neuroblasts is related to the prognosis of tumor patients, and EMT is closely related[28-30]. Sort1 is involved in the occurrence, development, and drug resistance of a variety of tumors, such as colon cancer, breast cancer and gastric cancer[31,32]. The expression of TCHH was correlated with immune infiltration in Colon cancer and this result is based on TIMER database analysis[33]. Very little research has been done on the relationship between the genes included in the model and the infiltration of immune cells in tumor tissue. In addition, the risk score calculated by this model is positively correlated with PDCD2LG2 (PD-L1) and Cytotoxic T lymphocyte-4, suggesting that the model may be helpful to guide the treatment of immune examination sites.

Taha-Mehlitz et al[24] study showed that the mutation and overexpression of FGFR3 were related to the low immune score of bladder cancer. It has been reported that Treg cells are highly infiltrated into tumors, which may lead to immunosuppression of tumor microenvironment and increase tumor immune escape[9]. Based on our data and previous studies, we have obtained a hypothesis that the increase of FGFR3 mutation will lead to the increase of Treg infiltration (because the mutation frequency of FGFR3 is significantly increased in the high-risk group), which will lead to the immune escape of tumor cells and lead to the poor prognosis of patients. This hypothesis may be a good research idea. The prediction of drug susceptibility may be an unexpected finding of this study.

Based on our data analysis, this model related to Tregs and NK cells can predict the prognosis of BUC, correlate with clinical data, and predict the sensitivity of bladder patients to chemotherapeutic drugs. However, it is not clear whether the changes of these genes in tumor tissue are the cause or result of Tregs and NK cells participating in the formation of TME. Here we merely provide some reference data for future researchers.

We constructed a prognosis-related model for patients with bladder tumors. This model is based on the infiltration of Treg cells and NK cells in tumor tissue. In addition to judging the prognosis of bladder cancer patients, it can also predict how sensitive patients are to chemotherapy and immunotherapy. Furthermore, based on this model, patients were divided into high and low risk groups, and variations in genetic alterations were found between the high and low risk groups.

In summary, Treg cell and NK cell percolation in this model provides new targets for targeted therapy of BUC. This study does have certain restrictions, though. Because there are many types of bladder cancer, yet, other research aspects may be overlooked, due to the pursuit of "accuracy" in the study of differentially expressed genes. In follow-up studies, we would consider looking at other subtypes to determine if there are any changes and investigating other immune-related processes for prognosis-related factors.

The role and function of Tregs and NK cells need to be confirmed using more sophisticated methodologies and approaches as our analysis is based on the TCGA database, which is openly accessible but has been mined.

CONCLUSION

Prognosis-related models of bladder tumor patients, based on Treg and NK cell percolation in tumor tissue. In addition to judging the prognosis of patients with bladder cancer, it can also predict the sensitivity of patients to chemotherapy and immunotherapy. At the same time, patients were divided into high and low risk groups based on this model, and differences in genetic mutations were found between the high and low risk groups.

ARTICLE HIGHLIGHTS

Research background

Regulatory T cells (Tregs) and natural killer (NK) cells play an essential role in the development of bladder urothelial carcinoma (BUC).

Research motivation

To identify genes associated with bladder cancer prognosis.

Research objectives

To construct a prognosis-related model to judge the prognosis of patients with bladder cancer, meanwhile, predict the sensitivity of patients to chemotherapy and immunotherapy.

Research methods

In this study, we analyzed publicly available datasets from two databases (https://portal.gdc.cancer.gov/ and https://www.ncbi.nlm.nih.gov/geo/), both with authoritative data, to predict the diagnosis and prognosis of the disease by bioinformatic analysis.

Research results

The stage and risk scores are independent prognostic factors in patients with BUC. Mutations in FGFR3 lead to an increase in Tregs percolation and affect the prognosis of the tumor, and additionally, EMP1, TCHH and CNTNAP3B in the model are mainly positively correlated with the expression of immune checkpoints, while CMTM8, SORT1 and IQSEC1 are negatively correlated with immune checkpoints and the high-risk group had higher sensitivity to chemotherapy drugs.

Research conclusions

The model can predict the sensitivity of patients to chemotherapy and immunotherapy in addition to their prognosis for bladder cancers.

Research perspectives

Prognosis-related models of bladder tumor patients, based on Treg and NK cell percolation in tumor tissue. In addition to judging the prognosis of patients with bladder cancer, it can also predict the sensitivity of patients to chemotherapy and immunotherapy. At the same time, patients were divided into high and low risk groups based on this model, and differences in genetic mutations were found between the high and low risk groups.