We obtained publicly available RNA-seq data for bladder cancer patients and the patient clinical information and mutation data from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). We transformed the RNA-seq data (fragments per kilobase million, FPKM) values into transcripts per kilobase million (TPM) values before performing log normalization. We performed gene annotation using the Ensembl database (http://asia.ensembl.org/biomart/martview). The GSE48075 cohort, which contains 142 bladder cancer samples, was retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi). By removing the samples with a vital state listed as “NA” using the “filter” function in R version 4.0.5 (https://www.r-project.org), we obtained 73 BLCA samples. We obtained a panel of gene signature matrices that distinguished 22 immune cell phenotypes, referred to as LM22, from the available literature (Newman et al., 2015). We also collected a set of 50 hallmark signatures (Liberzon et al., 2015) and a panel of biological function signatures deposited in the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) to conduct a GSVA analysis.

We analyzed the PD-1 and PD-L1 mutational signatures in the transcriptome and mutation profiles. For each gene out of a panel containing 16,542 genes, we regrouped 408 patients into mutant and wild-type cohorts. We carried out the Wilcoxon test to observe the differences in PD-1/PD-L1 expression between the mutant and wild-type groups and computed the expression fold change values between these two groups. The genes were only considered PD-1/PD-L1 mutational signatures under the following conditions: 1. mutational frequency >2%; 2. p value of Wilcoxon test <0.05; and 3. FC >1.5.

We classified the 408 bladder cancer patients into two groups (high and low expression) using RNA-seq count data. The top 25% PD-1/PD-L1 expression values were treated as a cut-off to differentiate the “high” and “low” expression cohorts. We installed a DESeq2 R Bioconductor package to analyse differentially expressed genes (DEGs) between the two cohorts. The significant DEGs met the criteria of |log2FC|>1.5 and an adjusted p value < 0.05. We used the “Enhanced Volcano” R package to draw volcano plots.

We combined PD-1/PD-L1 mutational signatures with PD-1/PD-L1-associated DEGs to analyze PD-1/PD-L1 genomic signatures using the “draw_venn” function in R.

Based on the 401 TCGA BLCA samples with available RNA-seq data and clinical information, we performed a Kaplan-Meier analysis of the overall survival (OS) rates to screen genes related to prognosis from the PD-1/PD-L1 genomic signatures. Patients were stratified by gene expression values greater than or below the median score. We applied the log-rank test to observe the differences between the two groups using the “survdiff” function in the “survival” R package, with p < 0.05 as the significance threshold. Subsequently, the GSE48075 population was used to test these selected genes for prognostic significance. We adopted the “survminer” R package to plot survival curves.

We analyzed the relationship between the SORL1 expression level and BLCA phenotypes (pathologic T, pathologic N, pathologic M, pathologic stage, primary therapy outcome success, sample type) using the online platform UCSC Xena (http://xena.ucsc.edu/).

We determined the protein expression of SORL1, both in BLCA and normal tissues, using IHC retrieved from the Human Protein Atlas (HPA, https://www.proteinatlas.org/) (Uhlen et al., 2010).

We analysed the gene expression correlation using the online platform cBioPortal (http://www.cbioportal.org/), which comprises various cancer genomic data. We analysed 408 mRNA expression datasets derived from bladder cancer samples (TCGA, Cell 2017). Additionally, a scatter plot was used to visualize the correlation between the log2-normalized gene expression value of SORL1 and PD-1/PD-L1, and the p value was used to determine statistical significance.

We used the Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) method (Newman et al., 2015; Zhang Z. et al., 2021) to quantify the infiltration level of 22 different immune cells. Contrasting the BLCA RNA-seq FPKM matrix and LM22, we estimated the abundances of immune cell types in 401 BLCA patients using the CIBERSORT R script obtained from the CIBERSORT platform (https://cibersort.stanford.edu/) after 1,000 permutations. We excluded the samples with a p value for the deconvolution algorithm greater than or equal to 0.05 and retained the remaining 204 samples for further analysis. We conducted a Wilcoxon test to evaluate the differences in immune infiltration levels between the SORL1 high and low expression groups.

We applied the Gene Set Variation Analysis (GSVA) R Bioconductor package to the reference gene sets to estimate the pathway enrichment scores in different samples (Hänzelmann et al., 2013). We used the limma R package to perform a differential analysis between high and low expression groups divided by the median SORL1 expression level. We defined the significantly differentially enriched signature based on adjusted p-values < 0.001 and |log FC|>0.3.

Then, we performed a Kaplan-Meier analysis and a log-rank test to evaluate the association between the selected differentially enriched signatures and BLCA survival outcomes. p < 0.01 was considered statistically significant.

We obtained an epithelial-mesenchymal transition (EMT) signature enrichment score (defined as the EMT score) from the previous step and computed the M2 macrophage immune infiltration abundance using the CIBERSORT method. We generated the scatter plot of the correlation between the EMT scores and M2 infiltration abundance through the “ggplot2” function in R. We estimated the linear regression line using the “lm” function and estimated the Spearman correlation using the “stat_cor” function in R.

To identify immune-relevant biomarkers that can predict the prognosis of bladder carcinoma, we explored PD-1/PD-L1-relevant genomic profiles in bladder cancer. Figure 1 shows the flowchart of this study.

We selected 408 bladder tumour patients out of 433 TCGA BLCA samples. These 408 samples were accompanied by both transcriptomic and mutation profile data. We then evaluated the mutation state of 16,542 genes in the selected samples. Their mutation state was studied with a focus on PD-1- and PD-L1-related genomic profiles. For each gene, we regrouped samples into mutant and wild-type cohorts. With the criteria that the mutant group to wild-type ratio had a mean PD-1/PD-L1 expression value >1.5 and a mutational frequency >2%, we obtained the PD-1/PD-L1 mutational signatures listed in Figures 2A,B.

Next, using a minimum |log2FC| of 1.5 and adjusted p value of no more than 0.05, we performed a differential expression analysis between the PD-1/PD-L1 high-expression group (102 samples) and low-expression group (306 samples). As a result, we obtained a panel of DEGs out of 53,170 genes, which are described in Figures 2C,D.

After Venn analysis, we obtained the PD-1 genomic signatures (n = 8) by calculating the intersection of the PD-1 mutational signature and the PD-1-related DEGs (Figure 3A). We carried out the same method for PD-L1 and obtained the PD-L1 genomic signatures (n = 14) (Figure 3B). Finally, we combined (n = 20) the PD-1- and PD-L1-related signatures (Figure 3C and Table 1).

Based on the clinical information obtained from the TCGA BLCA samples (401 cases), we performed the Kaplan-Meier survival analysis on 20 genes. Fortunately, we found that the genes TOX3 (p = 0.0043) and SORL1 (p = 0.042) were associated with BLCA (Figures 4A,B). To verify our finding, we applied Kaplan-Meier analysis to the GSE48075 cohort and found that the prognostic value of TOX3 lacked statistical significance (p = 0.17), but that BLCA patients in the SORL1 high-expression group showed significantly (p = 0.0065) more favourable OS outcomes than those in the SORL1 low-expression group in the GEO validation cohort (Figures 4C,D). In summary, SORL1 was identified as the only good prognostic indicator for BLCA overall survival outcomes.

After identifying the survival-associated gene SORL1, we verified the relationship between the expression levels of SORL1 and the clinical parameters. From Figure 5A, we found there was a rough trend that the lower the level of SORL1 expression, the later the T staging (p < 0.01). However, there was no statistical significance in N staging (p > 0.05) (Figure 5B). Figure 5C shows the association between the SORL1 expression level and the pathological M stage of bladder carcinoma. A comparatively lower SORL1 expression level indicated a greater risk for metastasis (p < 0.05). Figure 5D depicts a higher expression level of SORL1 in the earlier stages of bladder cancer (p < 0.05). Figure 5E demonstrates the influence of the SORL1 expression level on primary therapeutic outcomes. However, the influence did not achieve statistical significance (p > 0.05). Figure 5F shows that the mRNA expression level of SORL1 was significantly reduced in tumours (p < 0.01). Furthermore, the protein levels of the SORL1 gene were lower in tumour tissues (Figure 5G) than in normal tissues (Figure 5H). Overall, these findings confirmed that the downregulation of SORL1 expression is associated with worse survival outcomes in BLCA patients.

To investigate the association between the expression of SORL1 and the immune state, we performed a correlation analysis between SORL1 and PD-1/PD-L1 expression. We observed a negative correlation between the PD-1 (also named PDCD1) expression value and the PD-L1 (also named CD274) expression value with SORL1 expression. Additionally, the coefficients of the Spearman correlation were −0.32 (p = 2.83e−11) and −0.49 (p = 5.00e−26) (Figures 6A,B).

To clarify the survival mechanism related to the relationship between BLCA and the expression of SORL1, we explored the differences in tumour-infiltrating immune cells between patients with high and low SORL1 expression (the group allocation of patients was the same as the previous step). The results are depicted in Figures 7A,B. Figure 7C shows the landscapes of the proportions of 22 immune cells in the selected BLCA patients. We found that the differences in naïve B cells (p < 0.001), plasma cells (p < 0.001), activated memory CD4 T cells (p < 0.001), regulatory T cells (p < 0.001), resting NK cells (p < 0.001) and M2 macrophages (p = 0.002) were statistically significant. Additionally, the differences in naïve B cells, activated memory CD4 T cells, and M2 macrophages were much greater than the others (Figure 7B).

To investigate the SORL1-related pathway signatures, we performed GSVA on the TCGA bladder cancer samples.

Using GSVA, we obtained the enrichment scores of 50 hallmark signatures approximately in a normal distribution. We presented the association between the hallmark signatures and SORL1 expression level in the form of a heat map (Figure 8A). We identified three significantly differentially enriched signatures between the SORL1 high-expression and low-expression groups using the criteria of adjusted p-values < 0.001 and |log FC|>0.3. These were HALLMARK_ALLOGRAFT_REJECTION, HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION and HALLMARK_INFLAMMATORY_RESPONSE. Subsequently, the Kaplan-Meier analysis showed that the pathway signature EMT was associated with significantly worse survival outcomes (p = 0.0077) (Figure 8B).

For a set of biological process signatures, we also obtained 251 significantly differentially enriched signatures using the same criterion as the one described above. Figure 8C shows the top 50 differentially enriched pathways. We also applied Kaplan-Meier analysis to these 251 signatures and found that 14 signatures were associated with BLCA patient prognosis (Table 2). When we focused on the results, the signature regulation of macrophage fusion caught our attention. This signature was associated with poor prognosis among bladder cancer patients (p = 0.0067), as shown in Figure 8D.

To further explore the correlation between EMT and M2 macrophages in bladder carcinoma, we performed a correlation analysis on the EMT score and M2 infiltration abundance. Figure 8E depicts a significant positive correlation between M2 macrophage levels and the EMT signature (cor = 0.31, p = 1.2e−10).