Clinical neutrophil-related gene helps treat bladder urothelial carcinoma

Data sorting

In our study, the selected transcriptome data included the TCGA-BLCA and GSE31684 groups. The TCGA-BLCA cohort included 433 BLCA samples divided into 19 normal tissue sequencing samples and 414 tumor samples. After screening, 407 are left. We used somatic mutation data from TCGA database for internal validation. We, then selected the GSE31684 cohort from the Gene Expression Omnibus (GEO) database for external validation.

Landscape of infiltrating immune cells

Through CIBERSORT algorithm [13], we analyzed the sample data of TCGA-BLCA cohort patients to reveal the proportion of 22 tumor infiltrating immune cells (TIICs) in each sample.

Weighted gene co-expression network analysis

The purpose of using WGCNA is to help us build a gene relationship network to identify key genes [14]. We obtained the expression data of sample genes through the TCGA-BLCA cohort, and calculated the correlation coefficient between genes through the Person coefficient. Then, using the weighting function and adjacency function to obtain the parameters β to get the ideal soft threshold. The hierarchical clustering tree constructed by using TOM matrix is divided into 13 candidate modules with different colors by dynamic cutting algorithm, each module contains similar expressed genes. Finally, the relationship graph between 13 candidate modules and 22 kinds of immune infiltrating cells was obtained by correlation analysis between modules and phenotypic traits. Because this study is mainly aimed at neutrophils, we extracted the most relevant modules of neutrophils.

Cox and lasso regression analyses

Since the number of genes that affect prognosis screened out in this study is far greater than the study sample size, in order to avoid model distortion caused by multicollinearity among variables, Lasso+Cox survival analysis mode is used in this study [15]. To facilitate our exploration of the prognostic function of neutrophil-related genes, we screened 4,451 genes in the “grey” module. Firstly, we used univariate regression analysis to determine the genes related to survival. Then, we use lasso regression analysis to avoid over-fitting. Finally, the results were analyzed by multivariate regression to screen out 10 neutrophil-related genes that affect the prognosis of BLCA patients.

Validation of the prognostic neutrophil-related signature

To obtain the RS, we calculated the gene expression associated with prognosis and regression analysis coefficients based on TCGA cohort data. The formula used is as follows:

The median RS value is used for grouping, namely HRG and LRG. Kaplan Meier curve was used to compare the prognosis difference between the two groups. Finally, the value of the prognostic coefficient was proved by the ROC curve.

Establishment and verification of the nomogram

Our nomogram survival plot combined RS and clinical variables to predict patient survival at 1, 3, and 5 years. The validity of the model was verified by using a calibration curve.

Gene set enrichment analysis

We used the c2.cp.kegg.v7.4. symbols collection to study the functional annotation using GSEA software.

Correlation between tumor mutation burden and risk score

Somatic mutation data were acquired from TCGA database. Waterfall charts of the HRG and LRG were drawn with the “maftools” package of R. The survival difference of BLCA patients was analyzed according to median TMB and RS.

Correlation of risk score and tumor immune microenvironment

In order to study the potential relationship between RS and TIICs, we utiliazed seven methods to evaluate the immune infiltration. The proportion of stromal cells and immune cells in BLCA was estimated by ESTIMATE algorithm, and the immune score and matrix score were calculated. Add the two scores to get the ESTIMATE score, which predicts tumor content as well as calculates tumor purity in each tumor sample. The correlation between RS and tumor immune microenvironment was studied by Spearman correlation.

Gene set variation analysis(GSVA)

Though the MSigDB database to conduct relevant pathway analysis [16]. We performed GSVA in the GSVA package for distributional pathway activity estimation, in order to facilitate the assessment of relative pathway activity in a single sample.

Prediction of immunotherapeutic effect

Immune checkpoints refer to programmed death receptors and their ligands, which are the key targets that inhibit the function of immune cells. The immunophenotyping score (IPS) was used to assess the efficacy of immune checkpoint inhibitor therapy and determine tumor immunogenicity. The IPS can calculate scores for four immune phenotypes, including antigen presentation, effector, inhibition, and checkpoint, and the combination of scores of the four immune phenotypes is the IPS z-score. The higher the score, the more immunogenic is the sample.

Prediction of chemotherapeutic effect

We established ridge regression model based on tumor cell line chemosensitivity genomics and TCGA gene expression profile to explore RS and chemosensitivity. The half-maximal inhibitory concentrations (IC50s) of four chemotherapeutic drugs (Imatinid, Cislpain, EHT.1864, and lenalidomide) in patients with BLCA were evaluated using a predictive algorithm.

Statistical analysis

The Wilcoxon test and Kruskal Wallis test [17] were used to compare two groups and more than two groups of variables. Survival analysis was performed using the Kaplan–Meier log-rank test. The Chi-square test was used for analysis of the association between RS and TMB, and Spearman analysis was used to calculate the correlation between coefficients. Two-sided p<0.05 was considered to indicate statistical significance. All statistical analyses were performed using R software (version 4.1.1).

Landscape of TIME in BLCA

We calculated the proportion of TIICs in each sample through CIBERSORT algorithm (Attachment 1: Table S1), revealing the landscape of tumor infiltration microenvironment in BLCA. Then, the relationship between TIME and normal tissues and tumors was discussed and described through comprehensive thermography (Figure 1A). In order to further reveal the potential relationship between these immune cells, we tested the correlation between 22 immune cells in the TIME landscape (Figure 1B). It is noteworthy that neutrophils are positively correlated with activated mast cells (r=0.20; p<0.05) and negatively correlated with regulatory T cells (r=−0.14; p<0.05).

Figure 1:

Subpopulation of 22 immune cell subtypes and (A) proportional heat mapof 22 TICs in each BLCA sample. (B) Intrinsic correlation of 22 infiltrating immunecell types in BLCA.

Establishment of the WGCNA network

We calculate the parameters β=9 through the weighting function and adjacency function in WCGNA, and there is a good linear relationship between the validation results of R² when the visible soft threshold is 9 in Supplementary Figure S1A. Then 13 gene modules with different colors were obtained through the TOM matrix algorithm (Supplementary Figure S1B). Each row in Supplementary Figure S1C has 22 TIC types, and each column has candidate modules of the trait carrier gene. Notably, among the 13 candidate modules, the grey module was associated with neutrophils (cor=−0.23, p=1e−06). Our focus in this study was on neutrophils; therefore, we extracted the genes in the grey module (Additional File 1: Table S2) for further studies.

Development of risk signature

Using TCGA-BLCA data, we obtained expression data and follow-up information to facilitate follow-up studies on the prognostic value of the candidate genes. Using univariate Cox analysis, 4,451 neutrophil-related genes were found to have notable prognostic values (p<0.05, Additional File 1: Table S3). In order to get over fitting results, we conducted Lasso regression analysis on 610 genes screened (Figure 2A and B). Then we further determined 10 neutrophil-related genes (TM4SF1, PCED1B, GJD3, ZNF474, APOL2, PXDNL, BAIAP2, BCL2L14, SIGLEC6, and FAM166B) through multivariate Cox regression analysis. And among these central genes, PCED1B, GJD3, APOL2, FAM166B and BCL2L14 are all useful prognostic indicators (all risk ratios [HRs] are <1, Table S4).

Figure 2:

(A) Variation curve of regression coefficient with log (λ) in lasso regression. (B) Ten-fold cross-validation for tuning parameter selection in lasso regression. Vertical lines are drawn from the best data according to the minimum criterion and 1 standard error criterion. (C) Kaplan–Meier curve analysis showing the difference in overall survival between high-risk and low-risk groups in TCGA-BLCA cohort. Kaplan–Meier curve analysis shows the difference in overall survival with respect to APOL2, (D) gene and BAIAP2, (E) gene between the high-expression and low-expression groups. (F) Univariate Cox regression results for overall survival. (G) Multivariate Cox regression results for overall survival.

We calculated the risk score of the 10 central genes in patients with BLCA. Calculation of RS was performed as follows:

Based on the median of RS, we divided the sample into high-risk score group (HRG) and low-risk score group (LRG), so as to facilitate the grouping identification of later research.

Validation of risk prognostic features

From the Kaplan-Meier survival curve, LRG is superior to HRG in survival results (Figure 2C). Taking APOL2 and BAIAP2 as examples, it can be seen in Figure 2D that the survival of APOL2 with high expression is better than that with low expression. In Figure 2E, it can be seen that BAIAP2 survives better under the condition of low expression. According to the results of the two charts, the samples were divided into two groups: high expression group (low risk group) and low expression group (high risk group). From the Kaplan-Meier survival curves in Figure 2D and E, it can be seen that APOL2 and BAIAP2 genes have positive reactions to the prognosis of BLCA patients, suggesting that there is a certain relationship between the screening of neutrophil-related genes and the prognosis of the disease, and its risk score is feasible to be used to construct the prognosis model. We carried out the same verification for the other hub genes (Supplementary Figure S2). Univariate Cox analysis showed that the HR of RS was 1.282 (95 % confidence interval [CI], 1.223−1.344; Figure 2F). Multivariate Cox regression analysis results (HR=1.24295 % CI=1.182−1.306; Figure 2G) showed the effectiveness of RS as an independent prognostic indicator of BLCA. These results show that these 10 hub genes have good clinical prediction ability.

The expression pattern, sample survival status distribution and risk score of 10 genes in TCGA-BLCA cohort are shown in Figure 3A–C. In addition, we used the samples in the GSE31684 queue for external verification (Figure 3D–F). These results all support the above point of view. In the two data samples, we can see that the low risk score has a better survival state.

Figure 3:

(A) Confirmation of prognostic risk scores in TCGA cohort. (B) Polygenic model risk score distribution in TCGA cohort. (C) Survival status and duration of patients with BLCA in TCGA cohort. (D) Confirmation of prognostic risk scores in the GEO cohort. (E) Polygenic model risk score distribution in the GEO cohort. (F) Survival status and duration of patients with BLCA in the GEO cohort.

Functional analysis of neutrophil-related genes

We identified the functional enrichment of hub gene expression based on GSEA. Through the Kyoto Encyclopedia of Genesand Genomes(KEGG) pathway enrichment analysis, we found that the main pathways involved in APOL2 are graft vs. host disease, antigen processing and presentation, and Leishmania infection (Figure 4A). We performed the same analysis for the other central genes. KEGG pathway enrichment analysis indicated that the overexpression of BAIAP2 interacted with the signal transduction pathway in small cell lung cancer (Figure 4B). High expression of BCL2L14 and GJD3 was related to the metabolism of xenobiotics through the cytochrome P450 signal transduction pathway (Figure 4C and D). High expression of FAM166B was associated with the systemic lupus erythematosus signal transduction pathway (Figure 4E). High expression of PCED1B, SIGLEC6, and TM4SF1 genes was related to the signal transduction pathway of the intestinal immune network for IgA production (Figure 4F), but high expression of SIGLEC6, and TM4SF1 resulted from their co-action with the signal transduction pathways of viral myocarditis and Leishmania infection, respectively (Figure 4G and H). Overexpression of PXDNL was related to the taste signal transduction pathway (Figure 4I), and high expression of ZNF474 resulted from the interaction between taste transduction and focal adhesion (Figure 4J).

Figure 4:

GSEA of samples with high and low expression of two hub genes. (A) Gene set of samples enriched in APOL2 expression collected using KEGG. (B) Gene set of samples enriched in BAIAP2 expression collected using KEGG. (C) Gene set of samples enriched in BCL2L14 expression collected using KEGG. (D) Gene set of samples enriched in GJD3 expression collected using KEGG. (E) Gene set of samples enriched in FAM166B expression collected using KEGG. (F) Gene set of samples enriched in PCED1B expression collected using KEGG. (G) Gene set of samples enriched in SIGLEC6 expression collected using KEGG. (H) Gene set of samples enriched in TM4SF1 expression collected using KEGG. (I) Gene set of samples enriched in PXDNL expression collected using KEGG. (J) Gene set of samples enriched in ZNF474 expression collected using KEGG.

Correlation of risk signature with clinicopathological variables

In Figure 5A, we can clearly observe the distribution difference of clinically relevant variables in HRG and LRG. From Figure 5B–G, we can see that LRG presents a more optimistic result than HRG in terms of disease progression and prognosis.

Figure 5:

Clinical significance of prognostic risk characteristics. (A) Heat map showing the distribution of clinical characteristics and corresponding risk scores in each sample. Incidence of subtypes of LRG/HRG clinical variables: (B) sex, (C) tumor grade, (D) clinical stage, (E) T stage, (F) N stage, and (G) M stage.

Construction of prognostic nomogram

The ROC curve shows that the area under the curve (AUC) values of 1 year overall survival rate (OS), 3 year OS and 5 year OS are 0.754, 0.776 and 0.796 respectively. These values are sufficient to explain the high prognostic effectiveness of the features analyzed in this study (Figure 6C). In addition, AUC was used to analyze 1 year (Figure 6D), 3 year (Figure 6E), and 5 year OS (Figure 6F) for a combination of risk characteristics, sex, and clinical stage. Finally, we drew a prognostic nomogram (Figure 6A) so that we could more intuitively evaluate the impact of risk characteristics, sex, TNM stage, and clinical stage on OS in patients with BLCA. The validation results from the calibration curve showed that the nomogram was predictive (Figure 6B).

Figure 6:

Evaluation of risk signal prediction efficiency. (A) Nomogram used to predict survival in patients with BLCA. (B) 1-, 3-, and 5-year survival nomogram calibration curves. (C) ROC analysis was used to assess the predictive value of prognostic features. The AUC for the risk score was used to predict overall survival at 1, 3, and 5 years and other clinical characteristics (D–F).

Association of risk signature with TMB

Next, we explored the potential associations between TMB and risk characteristics by comparing TMB in samples from the LRG and HRG, and we found that high TMB level is more closely related to LRG (p<0.001, Figure 7A). To facilitate the study of the effect of TMB on the OS of patients with BLCA, we constructed a Kaplan–Meier survival curve, and the grouping of samples from the high- and low-TMB groups in the curve was based on the median TMB value (p<0.001, Figure 7B). Subsequently, we analyzed the correlation between TMB and RS, as shown in Figure 7C, and we found that RS was inversely proportional to TMB. Figure 7D clearly illustrates the interaction between TMB and RS in the prognosis of patients with BLCA.

Figure 7:

Correlation between risk score and TMB. (A) Differences in TMB between the HRG and LRG. (B) Scatterplots depicting the positive correlation between risk scores and TMB. (C) Kaplan–Meier curves of the high- and low-TMB groups. (D) Kaplan–Meier curve stratification of patients according to TMB and risk signature. The oncoPrint was constructed using high risk score (E) and low risk score (F). TMB, tumor mutational burden; HRG, high-risk group; LRG, low-risk group.

Subsequently, based on the distribution and type of somatic gene mutations in different RS subpopulations, we created correlation maps of the respective somatic variants of the HRG and LRG (Figure 7E and F). The mutation profiles of significantly mutated genes revealed that TP53 (49 vs. 44 %) had more somatic mutations in the HRG core subtypes, whereas TTN (35 vs. 44 %), KDM6A (19 vs. 25 %), and FGFR3 (7 vs. 21 %) had more somatic mutations in the LRG.

Risk signature in TIME context

Spearman correlations were used to analyze the correlations between the risk signatures of neutrophil-related genes and TIME (Figure 8A, Table S5). According to the analysis results of ESTIMATE score, immune score is higher than stromal score in the LRG. However, the overall ESTIMATE score decreased (p<0.01, Figure 8B). After correlation validation using the four methods—CIBERSORT (Figure 8C), CIBERSORT−ABS (Figure 8D), QUANTISEQ (Figure 8E), and TIMER (Figure 8F)—our analysis was confirmed to be correct.

Figure 8:

Estimated abundance of tumor-infiltrating cells. The (A) high-risk group had a stronger correlation with tumor-infiltrating immune cells, as shown by the Spearman correlation analysis. (B) Association between prognostic risk signatures and central immune checkpoint genes. The correlations predicted using the four methods—CIBERSORT (C), CIBERSORT−ABS (D), QUANTISEQ (E), and TIMER (F)—were validated.

Enrichment of signaling pathways in the low/high-risk groups

Through GSVA analysis, we found that RS was negatively correlated with KEGG/PPAR signaling pathway, which indicated that the disease occurrence and development of LRG was related to this pathway. On the contrary, most signal pathways in Figure 9A are in positive proportion to RS.

Figure 9:

Enrichment pathways for GSVA. (A) Heat map showing the correlation of representative pathway items in KEGG with risk scores. Predicting immunotherapy response. (B) Association of expression levels of immune checkpoint blockade genes with risk scores. (C–F) IPS distribution map. Estimates of chemotherapeutic effect and risk scores. Sensitivity analysis of (G) EHT.1864 in patients with high and low risk scores. Sensitivity analysis of (H) lenalidomide in patients with high and low risk scores. (I) Sensitivity analysis of cisplatin in patients with high and low risk scores. (J) Sensitivity analysis of Imatini in patients with high and low risk scores.

Prediction of response to immunotherapy

Because of the lack of data on immunotherapy in TCGA-BLCA dataset, we could only use the expression levels of immune checkpoint blockade-related genes to evaluate the prognostic effect of immunotherapy. Moreover, we identified other immune checkpoint blockade-related genes (i.e., TNFRSF8, TNFRSF14, TMIGD2, and CD86) whose expression was significantly proportional to risk characteristics and determined some genes associated with immune checkpoint blockade, such as ICOS, CD27, CD244 and IDO2, whose expression weas inversely associated with risk-related traits (Figure 9B). It can be seen in Figure 9B that GJD3 is negatively correlated with most of the immune checkpoint block-related genes, while FAM166B is positively correlated with most of the immune checkpoint block-related genes, it seems contradictory by comparison, but this is not the case. The influence of neutrophil-related genes on risk score includes the occurrence and development of tumor and the effect of immunotherapy. We can from Figure 9C that the finding suggests that the LRG patients are more inclined to choose new immunotherapy with the ICIs and it can contribute to our understanding of the relationship between immunotherapeutic effects and RS.

Prediction of response to chemotherapy

After data analysis, IC50 value is obtained, which allows us to see the difference of therapeutic sensitivity of four chemotherapy drugs (Cisplatin, Imatini, EHT.1864, and lenalidomide) to HRG and LRG. The HRG was more sensitive to EHT.1864 (Figure 9G) and Lenalidomide (Figure 9H) than the LRG, while the LRG was more sensitive to Cisplatin (Figure 9I) and Imatini (Figure 9J). These results suggest that feasible strategies for treatment can be devised by exploring the relationship between risk characteristics and chemotherapeutic drug sensitivity.