A novel 6-gene signature derived from tumor-infiltrating T cells and neutrophils predicts survival of bladder urothelial carcinoma

Study design and datasets

As the flow chart shows (Figure 1), the samples of the TCGA-BLCA dataset were randomly separated into the training (n=197) and testing (n=196) sets for the construction and verification of prognostic model, respectively. The clinicopathological features of patients in the TCGA-BLCA dataset are listed in Table 1. Chi-squared test revealed that several potential confounding factors such as age, gender, tumor-node-metastasis (TNM) stage, tumor size, lymph node metastasis, and distant metastasis were all evenly distributed in the training and testing sets with P > 0.05. For the external validation of the identified model, a dataset combining three GEO datasets (GSE13507, GSE31684, and GSE48276) and a cohort containing 46 BUC patients recruited from JSPH were further analyzed. The clinicopathological features of patients in the two external validation sets are shown in Table 2. The study was approved by the Institutional Ethical Committee of JSPH (ID: 2016-SRFA-148) and written informed consent had been obtained from each participant.

Univariate Cox regression analysis in the TCGA-BLCA dataset revealed that the proportion of tumor-infiltrating T cells in TME was a favorable prognostic factor for BUC (hazard ratio [HR]: 0.975, 95% confidence interval [CI]: 0.963-0.988, P < 0.001) while neutrophils proportion was an unfavorable prognostic factor (HR: 1.082, 95% CI: 1.016-1.151, P = 0.013). As shown in Supplementary Table 2, a total of 164 genes whose expression levels were significantly correlated with both T cells and neutrophils proportions were defined as candidate genes for model construction (P < 0.05; spearman's rank correlation analysis).

Identification of a T cell and neutrophil-associated prognostic model

Among the 164 candidate genes, 10 genes [Epithelial membrane protein 1 (EMP1), Neuroblast differentiation-associated protein (AHNAK), Solute carrier family 1 member 6 (SLC1A6), Mast cell-expressed membrane protein 1 (MCEMP1), Ras guanyl-releasing protein 4 (RASGRP4), Family with sequence similarity 180 member A (FAM180A), Serine protease 8 (PRSS8), Six transmembrane epithelial antigen of prostate 4 (STEAP4), Heat shock protein family A member 1 like (HSPA1L), and ADAM metallopeptidase with thrombospondin type 1 motif 16 (ADAMTS16)] were remarkably correlated with the prognosis of BUC patients in the TCGA training set according to the results of univariate Cox regression analysis (P < 0.01; Supplementary Table 3). Lasso regression analysis was further conducted to reduce redundancy (Supplementary Figure 1). Stepwise multivariate Cox regression analysis further identified a 6-gene model with the maximum prognostic value which included EMP1, RASGRP4, HSPA1L, AHNAK, SLC1A6, and PRSS8 (Figure 2A). According to the constructed Cox proportional-hazards model, risk score of a BUC patient could be calculated by the formula: 0.194173005532207 * EMP1 expression + 0.300051365964018 * RASGRP4 expression - 0.25824698561501 * HSPA1L expression + 0.426348983855879 * AHNAK expression + 0.29902852359222 * SLC1A6 expression + 0.29565935174874 * PRSS8 expression. Multivariate Cox regression analysis revealed that except HSPA1L, the other 5 genes were all independent unfavorable prognostic factors for BUC patients (HR > 1; P < 0.01; Figure 2B).

In the TCGA training set, risk score of each patient was calculated. The median value of the risk scores (-0.186984284651394) was regarded as the cutoff value to distinguish high-risk and low-risk patients. Kaplan-Meier survival analysis revealed that in the training set, patients with higher risk scores had significant worse overall survival (OS) than those with lower risk scores (P < 0.001; Figure 3A). Generally, as the risk score increased, patients did have shorter survival time and higher mortality rate (Figure 3B). The multiple receiver operating characteristic (ROC) curves further indicated that compared with other clinical parameters, the estimated risk score [Area under the Curve (AUC) = 0.766] had higher prognostic value for the survival prediction of BUC patients (Figure 3C). In the TCGA testing set, the identified model proved to have stable prognostic performance for BUC. In this cohort, OS of the estimated high-risk patients was significantly worse than the estimated low-risk patients (Figure 3D, 3E). Moreover, risk score (AUC = 0.711) was still superior to other clinical factors for survival prediction in BUC (Figure 3F). Predictive value of the identified 6-gene model remained high when the data of the training and testing sets were combined (Supplementary Figure 2 and Table 3). In addition, expression levels of the 6 genes in tumor tissues became generally higher as risk scores increased.

Construction of a comprehensive nomogram

Based on the TCGA-BLCA dataset, a nomogram combining risk score and other clinical parameters including age, gender, and TNM stage was constructed to facilitate survival probability prediction for BUC patients (Figure 4A). Calibration curve analysis confirmed that the theoretical 3-year and 5-year cancer specific survival rates were in good agreement with the actual survival rates (Figure 4B, 4C), which demonstrated the high effectiveness of the identified signature.

Clinical- and immune-correlation analysis

Subgroup analyses and Spearman's rank correlation analyses were performed in the TCGA-BLCA dataset to explore the association between the identified signature and clinical parameters as well as tumor immunity. As shown in Figure 5, higher risk score was associated with the condition of patient’s elder age (Figure 5A), later TNM stage (Figure 5B), higher invasion depth (Figure 5C), and more lymph node metastasis (Figure 5D). For each tumor sample, the immune score and stromal score were calculated respectively using the ESTIMATE algorithm to quantify the proportions of immune and stromal cells in TME. Estimated risk scores were in positive correlation with the contents of both immune cells (r = 0.298, P < 0.001) and stromal cells (r = 0.416, P < 0.001) in BUC tumor tissues (Figure 5E, 5F). Moreover, higher risk score was significantly correlated with less T cells (r = -0.321, P < 0.001) and more neutrophils (r = 0.345, P < 0.001) infiltration in TME (Figure 5G, 5H).

Expression levels of the 6 key genes were also compared between tumor tissues and adjacent normal tissues in TCGA-BLCA dataset. The mRNA expression levels of AHNAK, EMP1, RASGRP4, and HSPA1L were significantly down-regulated in BUC tumor tissues compared with normal tissues while SLC1A6 and PRSS8 were significantly up-regulated (P < 0.05; Supplementary Figure 3). Moreover, gene expression levels were further compared between BUC patients aged < 70 and ≥ 70 years. Only HSPA1L showed significantly lower expression in tumor tissues in elder BUC patients compared with the younger counterparts (P = 0.021; Supplementary Figure 4); there were no significant expressing differences for the other 5 genes. The typical immunostaining graphs of these genes in BUC tumor tissues were provided by the human protein atlas (HPA) database (Supplementary Figure 5).

Functional exploration by gene set enrichment analysis (GSEA)

In the TCGA-BLCA dataset, the differential pathways and biological processes enriched in the high-risk and low-risk groups were listed in Supplementary Table 4.1, 4.2 (false discovery rate [FDR] < 0.05). A number of immune-related pathways (i.e., ‘cytokine secretion’, ‘interleukin-6 production’, ‘leukocyte migration’, ‘cytokine-cytokine receptor interaction’) were identified. Figure 6 lists some of the significant immune-related pathways involved in the high-risk groups compared with the low-risk groups. The results further demonstrated the close relationship between the identified signature and the immune status of BUC.

External validation in GEO datasets

In the combined GSE13507, GSE31684, and GSE48276 dataset, risk scores of the 331 BUC patients were calculated using the same formula, and patients were divided into high- and low-risk groups based on the established cutoff value. OS of the predicted high-risk group was significant worse than that of the low-risk group (Figure 7A, 7B). Consistent with the TCGA-BLCA dataset, patients’ risk score had significantly negative correlation with T cells proportions in TME; while for neutrophils, the relationship was reversely positive (Figure 7C, 7D). It was further confirmed that tumor-infiltrating neutrophils and T cells might be unfavorable and favorable prognostic factors, respectively, on the long-term survival of BUC patients. For BUC, the identified T cell and neutrophil-associated 6-gene signature had potent prognostic value.

Prognostic role of the 6 key genes in JSPH cohort

Expression levels of the 6 genes in 46 BUC tumor tissues were analyzed by quantitative reverse transcription polymerase chain reaction (qRT-PCR) in the JSPH cohort. After z-score normalization of gene expression levels, risk scores of BUC patients in the JSPH cohort were calculated according to the identified model. As shown in Figure 8A, in the JSPH cohort, patients with higher risk scores had significant worse OS than low-risk patients (P < 0.05; Kaplan-Meier survival analysis). For each of the 6 genes, BUC patients were divided into high- and low-expression groups according to median expression values, and the OS rates of the two subgroups were further compared by Kaplan-Meier survival analysis. As shown in Figure 8B, patients with higher AHANK, EMP1, SLC1A6, and RASGRP4 expression had significant worse OS than the corresponding low-expression group with P < 0.05.

Data source and processing methods

The survival time and clinicopathologic information of all the patients with primary BUC was downloaded from TCGA-BLCA project. Subjects with follow-up period less than 30 days were excluded. All the 393 samples in the TCGA-BLCA dataset were randomly divided into the training and testing sets following the principle that clinicopathologic features including age, gender, and TNM stage were evenly distributed in the two groups to avoid potential confounding effects. Tumor transcriptome profiling data including raw counts and Fragments Per Kilobase of transcript per Million (FPKM)-normalized RNA-seq data were obtained using the gdc-client tool (https://portal.gdc.cancer.gov/). The human reference genome assembly GRCh38 was used to annotate protein-encoding genes. The transcriptome expression profiles of 18,321 protein-coding genes were normalized using ‘limma’ R package [35] and integrated into a matrix for future analysis.

GSE13507, GSE31684, and GSE48276 datasets with complete gene expression information and overall survival time of BUC patients were downloaded from GEO (https://www.ncbi.nlm.nih.gov/) for external validation. The corresponding platform files (GPL6102, GPL570, and GPL14951) of each array were further downloaded for the annotation of gene symbols. In total, 331 BUC patients from the three datasets with follow-up time more than 1 month were combined to enlarge the sample size. The ‘sva’ [36] and ‘limma’ [35] R packages were used to remove batch effect and normalize expressing data, respectively. The gene expression profiles and survival information of the three datasets were thus combined into one matrix for further analysis. The z-score normalization method was used to compare the different data types of RNA-seq and microarray.

Proportion calculation of tumor-infiltrating immune cells

The ‘CIBERSORT’ algorithm was applied to calculate the proportions of tumor-infiltrating leukocytes in BUC tumor tissues [10]. According to the developers’ instruction, the ‘CIBERSORT’ analytical tool (https://cibersort.stanford.edu/) could accurately quantify cell fractions by the deconvolution of bulk tumor gene expression profiles rather than depend on traditional cell enumeration methods such as flow cytometry and immunohistochemistry. The signature file ‘LM22’ consisting of characteristic genes that specifically distinguish 22 mature immune cells was downloaded for reference. The normalized gene expression matrix of the TCGA-BLCA dataset and the combined GEO dataset was separately uploaded to the online application to generate the corresponding leukocyte proportion matrix. For each sample, a deconvolution P-value was calculated to quantify deconvolution confidence. Only samples with P-value < 0.05 were considered for further analysis. In addition, the fraction of immune and stromal cells in tumor tissues was further estimated using the ‘ESTIMATE’ algorithm following the author’s instruction [11]. The calculated ‘Immune score’ and ‘Stromal score’ represented the relative proportions of tumor-infiltrating immune and stromal cells, respectively.

Selection of candidate genes

The proportion of overall tumor-infiltrating T cells in each BUC sample was calculated by adding together the percentages of ‘CD8⁺ T cells’, ‘naïve CD4⁺ T cells’, ‘resting memory CD4⁺ T cells’, ‘activated memory CD4⁺ T cells’, ‘follicular helper T cells’, ‘regulatory T cells (Tregs)’, and ‘gamma delta T cells’ according to the results of CIBERSORT deconvolution. Univariate Cox regression analysis indicated that the proportions of T cells and neutrophils in BUC tumor tissues were significantly associated with the OS of patients in the TCGA-BLCA dataset (P < 0.05), which gave rise to the possibility of identifying a T cell and neutrophil-associated prognostic signature for BUC. Therefore, candidate genes for model construction were screened out following the three main criteria: (i) gene expression levels were significantly correlated with the proportions of T cells and neutrophils in tumor tissues (P < 0.05, spearman's rank correlation analysis); (ii) the average read count of the gene was more than 10 to ensure adequate abundance for detection; (iii) the gene was significantly associated with the prognosis of BUC patients (P < 0.01, univariate Cox regression analysis).

Construction of a cox proportional-hazards model and nomogram

The TCGA-BLCA training set was analyzed for the construction of prognostic model. In order to minimize overfitting, Lasso regression analysis was firstly performed for variable selection using the ‘glmnet’ R package [37]. Stepwise multivariable regression analysis was then conducted to determine the covariates in the Cox proportional-hazards model. The risk score of each patient in the training, testing, and external validation sets was calculated separately according to the coefficients estimated by the model. The median of the risk scores in the training set was regarded as the cutoff value to determine high-risk and low-risk patients. In addition, a nomogram combining the estimated risk score and several clinical parameters was further constructed in the TCGA-BLCA dataset to facilitate prognostic prediction of BUC patients.

Functional enrichment analysis

GSEA was performed to explore the underlying function of the identified signature. The transcriptome profiling data was compared between high-risk and low-risk patients in the TCGA-BLCA dataset. The classical gene sets of Kyoto Encyclopedia of Genes and Genomes pathways and Gene Ontology project were analyzed to decipher the phenotypic differences between the two groups. For each analytical process, enrichment score (ES) and significance of ES were estimated, and normalized enrichment score and FDR were further calculated to examine the results of functional enrichment analyses. An FDR cutoff value of 5% was considered in this test.

QRT-PCR verification in tissue samples

Paraffin-embedded tissue specimens were collected from a total of 46 BUC patients who underwent operations. With the approval of the Institutional Ethical Committee (ID: 2016-SRFA-148), these patients were recruited from JSPH during 2011 and 2017 according to agreement principle. Patients were actively followed for survival information. Total RNA was extracted from tissue specimens using RecoverAll™ Total Nucleic Acid Isolation Kit (Ambion, Austin, TX, USA). Reverse transcription and qRT-PCR reactions were performed using PrimeScript™ RT reagent Kit (Takara, Kyoto, Japan) and SYBR^® Premix Ex Taq II (Tli RNaseH Plus) (Takara, Kyoto, Japan) following the manufacturer’s protocols, respectively. The sequences of PCR premiers used in this study were listed in Supplementary Table 1. GAPDH was considered as a reference gene for normalization, and the 2^-∆∆Ct method was used to analyze the relative expression of target genes [38]. Kaplan-Meier survival analysis was conducted to assess the association between gene expression levels and BUC prognosis.

Statistical methods

Chi-squared test was performed to determine whether the distribution of clinicopathologic characteristics was balanced between the TCGA training and testing sets. The clinical outcomes between patients with high and low risk scores were compared by Kaplan-Meier survival analysis. ROC curve analyses were conducted to assess the discriminating capability of the identified markers and make comparison with other clinicopathologic parameters. Calibration curve analysis was further considered to measure the predictive performance of the model. Multivariate Cox regression analysis was conducted to identify independent prognostic factors. Risk scores were compared between different subgroups using Mann-Whitney U test or Kruskal-Wallis H test. The correlation between risk score and immune cell fraction was assessed by Spearman's rank correlation analysis. R3.6.5 (https://www.r-project.org) and SPSS 25.0 (SPSS Inc., Chicago, IL, USA) software were applied for data analysis and graphing. A two-tailed P-value less than 0.05 was statistically significant.

Data availability

The data used to support the findings of this study are included within the article.