Gene expression and coexpression alterations marking evolution of bladder cancer

Dataset mining

A comprehensive data mining strategy was employed to retrieve studies applying -omics technologies in BLCA. The overall workflow is summarized in Figure 1. All genomic urothelial cancer data from cBioportal (including The Cancer Genome Atlas) were downloaded (5/1/2020). Gene Expression Omnibus (GEO) was queried for transcriptomics, additional genomics or protein array datasets using the search terms “bladder cancer” and “urothelial carcinoma”. We also queried ArrayExpress using the special filter “Array express data only” to obtain any additional datasets missing from GEO. All cohort data published or updated between 2010 and 2019, annotated as Homo sapiens, coming from tissue samples with sample size ≥10, were initially retrieved (25/1/2020). All used datasets were published and downloaded anonymized.

• Download figure
• Open in new tab

Figure 1: Study design and workflow for the integrative analysis of primary BLCA transcriptomes.

Integration of the transcriptomes

Microarray data were summarized to gene level with the package oligo [18]. Affymetrix were normalized with the RMA method, while Illumina were filtered based on detection p value < 0.05, followed by quantile normalization, addition of 1, and transformation to the natural logarithmic scale. Data were annotated using biomaRt (v2.42.1) [19] and mircroarray probes matching to multiple genes were excluded from downstream analysis. The probe with the highest mean across arrays was selected as representative in cases where multiple probes matched to the same gene. Merging of expression matrices was based on the Hugo Gene Symbol using the intersection of genes between studies.

Adjustment for batch effects

To correct for batch effects across different studies in the discovery data, ComBat [20], removeBatchEffect [21], and naiveRandRuv [22] were evaluated. ComBat performed best and was chosen for further use (manuscript in preparation). The quality of corrected data was assessed with BatchQC [23] (Supplementary Figures 1-2), with boxplots of expression distribution per sample and principal component analysis plot of sample relationships, with gene expression comparisons of housekeeping to other genes, and with a set of 12 BLCA markers with known regulation across normal-NMI-MI or across normal-low grade-high grade disease (Supplementary Figure 3).

Differential expression, monotonicity and functional annotation

To identify genes that form a continuum of changes across BLCA stages, each disease stage was compared against normal samples. Genes significantly affected (Mann-Whitney p < 0.05 and had also same orientation of change in all comparisons), were extracted (n = 3,108 genes). We refer to this set as Concordantly Differentially Expressed Genes (CDEGs). Monotonicity for a gene was defined as being a CDEG and additionally having a continuously larger/smaller fold change with increasing stage (fold change, as defined in each disease stage versus normal control). Functional annotation and enrichment were performed using PubMed and the online tool GeneCards (https://ga.genecards.org/), respectively.

Monotonicity in pathway activation and analysis for gene coexpression

We utilized CDEGs to infer pathway activation scores and to create stage coexpression networks. Pathway activation scores per sample were calculated with the ssGSEA-GSVA method [24], using the Molecular Signature Database libraries of Hallmark, Canonical Pathways (Reactome subset), C3 (GO biological processes subset) and C5 (GTRD subset of transcription factor targets). Dorothea (https://github.com/saezlab/dorothea) [25] was utilized to assess regulon activity. Activation scores across stages were compared with Mann-Whitney tests, and direction of change was defined based on fold change (= Mean of stage – Mean of T0). Pathways showing a monotonal change in the activation scores with increasing stage were extracted. Monotonicity for a pathway was defined similarly to the previous section, i.e. being significantly different in all stage comparisons to T0, and also having a continuously larger/smaller fold change with increasing stage. Stromal infiltration scores were imputed with the ESTIMATE algorithm [26].

Gene-pair coexpression weights (non-negative) were approximated with ensemble learning, using GENIE3 [27] and direction of coexpression (positive/negative) was determined by Spearman’s coefficient. Stage specific networks were constructed with igraph using the top 5,600 positively correlated gene pair interactions with the highest weights for each of the stages. Networks were analysed with Louvain clustering [28] to identify local modules of significantly higher coexpression relatioships (communities) within each stage. The top 5 in size (= number of genes) communities of coxpressed genes per stage were analysed for Gene Ontology Biological Processes [29]. Significance was thresholded at p<0.05 after FDR adjustment according to Benjamini-Hochberg. Hub genes in the monotonal networks were defined based on the betweenness centrality metric.

Statistical analysis

For continuous variables, significance was defined at Mann-Whitney p < 0.05, unless stated otherwise. Categorical variables were investigated for significance with the Chi-squared test and were adjusted for multiple hypothesis (CRAN package RVAideMemoire). All statistical tests were two sided. Gene Ontology and Reactome pathway enrichment analysis was conducted with functions from the package ClusterProfiler that performs a two-sided hypergeometric test. All reported correlation scores correspond to the Spearman’s Rank coefficient. Kaplan-Meier analysis and cox proportional hazards regression were performed with the CRAN packages survminer and survival, and statistical significance was determined with the log-rank method. CIBERSORT analysis was conducted in the web platform https://cibersort.stanford.edu/, and only samples with successful deconvolution (p < 0.05, n = 350 samples) were further used for the comparisons of relative immune populations among stages. All processing, analyses and visualizations were created in the R programming language (version 4.0.2), using predominantly Bioconductor libraries.

Cohort description and data quality controls

To perform a comprehensive investigation of publicly available omics studies, we collected 105 datasets comprising more than 8,000 individuals (Figure 1). Minimum inclusion criteria for datasets were defined as having at least stage information per subject. To maintain high integrity, we selected microarray data quantified by commercially available single-color channel vendors (Affymetrix and Illumina). Samples or datasets collected after administration of neoadjuvant chemotherapy (NAC) or not clearly annotated as primary were excluded, generating a final dataset corresponding to 1,508 patients. This included in total 12 microarray studies which were employed as a discovery set (Table 1) and in addition, the TCGA 2017 BLCA project based on RNA-seq analysis was used as a validation set (Figure 1).

The compiled discovery cohort comprised 1,058 primary bladder cancer tumor transcriptomes of treatment-naïve patients without prior cancer history, along with 81 adjacent to the cancer bladder samples, a total of 1,139 gene expression profiles (Table 1). The ratio of men:woman was 3.5:1, with equal distribution among NMI and MI disease (p = 0.99) and similar mean ages at baseline diagnosis (68 years, p = 0.81). Percentages of normal samples, NMI, and MI in the dataset were 7.1%, 43.5%, and 49.4%, respectively (Table 1), with the grade distribution being as follows: 16.5% low grade, 48.8% high grade disease with the remaining samples lacking available grade information. Detailed histological records were missing for 71.5% of the cohort, with the most frequently reported histology among the available records being urothelial/papillary (= 23.3%) with squamous differentiation being the most frequent variant (= 1.3%). Annotations included mutation status for FGFR3, PI3KCA, RAS, RB1 and TP53.

Removing batch effects with ComBat was assessed with Relative Log Expression and Principal Component Analysis plots (Figure 2A), as well as with comparative analysis of expression levels between houekeeping and non-housekeeping genes (Figure 2B), with the algorithm BatchQC (Supplementary Figures 1 and 2) [23], or using a set of 12 positive BLCA markers with known regulation across normal-NMI-MI or normal-low-high grade disease (Supplementary Figure 3). ComBat successfully eliminated batch effects while maintaining the biological variability (Figures 2A-C). BatchQC reports indicated that the variability in the corrected data was generally explained by stage rather than the batch variable (Supplementary Figure 2), allowing for an in-depth downsteam analysis.

• Download figure
• Open in new tab

Figure 2:

ComBat corrected expression data of the discovery cohort and respective differential gene expression analysis A) Relative Log Expression (right) and Principal Component Analysis plots (left) showing gene expression levels per sample and sample relationships, respectively; the different colors denote the different datasets B) dispersion-to-expression plots illustrating preservation of housekeeping gene properties (higher mean expression and lower median absolute deviation (MAD) compared to non-housekeeping genes) in the corrected data, regardless of microarray (top) or RNAseq (bottom) defined housekeeping genes; C) numbers of differentially expressed genes (DEGs) between disease stages and normal adjacent urothelium (top), along with respective group sizes (bottom), E) heatmap of the top (based on fold change) 20 downregulated and top 20 upregulated monotonically changing genes with increasing disease stage.

Monotonically changing genes in BLCA regulate growth and cell-cycle progression

For initial assessment of the gene expression relations between normal adjacent urothelium (referred to as T0) and cancerous samples, we performed differential expression analysis of T0 versus NMI and T0 versus MI cancers. When compared to T0, a strong positive correlation (r=0.83, p < 2.1×10^-16) was observed between the two sets of log fold change values, suggesting mutual concordance in abundance change and directionality between NMI and MI tumors. To investigate transcriptional changes associated with increasing malignancy, CDEGs were extracted from the comparisons Ta vs T0, T1 vs T0, T2 vs T0, T3 vs T0 and T4 vs T0. A total of 157 genes were identified having a monotonal (i.e. continuously increasing or decreasing in comparison to T0) change in expression across stages, of which 118 were up- and 39 were downregulated with increasing stage (Table 2). Functional annotation revealed that for 46 of these genes, experimental evidence on mediating cell cycle progression exists (Table 2; Supplementary Table 1). Upregulated cell-cycle associated genes (n = 44) were not phase specific and included cyclins, DNA polymerases, regulators of the cohesin complex and kinetochore components. The two downregulated cell-cycle associated genes, BTG2 and CDC14B, are both tumor suppressors. The list included 23 genes involved in signal transduction (Table 2), 6 of which (ARHGAP11A, AURKA, CDKN3, PBK, PLK1, RRM2), promote cell-cycle progession and were all upregulated with increased stage. The data also indicated an overactivation of the Wnt pathway with increasing disease stage, with its upstream inhibitor APCDD1 being downregulated and its activating ligand WNT2 upregulated (Supplementarty Table 1). Fourteen of the 157 genes were transcriptional or translational regulators (Table 2), including genes with known upregulation in bladder cancer (transcription factors E2F1, DEPDC1 [30, 31]). Based on the monotonal changes with increasing stage, increased androgen receptor activity may be predicted, as both its translational enhancer BUD31 [32] and its downstream transcritption factor ELK1 [33] were upregulated. Four (HTR2C, LRP8, NENF, NMU) of the 157 genes are involved in neurotransmission or neuronal development, all upregulated (Table 2). Among the 157 genes, 21 genes were of not well described or unknown function (Supplementary Table 1), including the oncogenic factor TRIM65 [34] found upregulated with bladder cancer stage. Additionally, upregulation was detected for the cisplatin resistance gene CLPTM1L, as well as for SRD5A1, which catalyzes the conversion of testosterone and progesterone to their 5-alpha forms [35]. Further functional enrichment using GeneCards for the 157 genes verified their involvement in cell-cycle pathways, with top hits being related to the regulation of the Anaphase-promoting (APC) complex (score = 31.53), to PLK1 (score = 24.47) and Aurora B (score = 20. 95) signaling, as well as to TP53 (score = 19.06) and RB1 (score = 17.85) cell cycle checkpoint control (Supplementary Table 2). Using the same tool, analysis for compounds with potentially therapeutic benefit provided as top hit the DNA alkylating agent Bendamustine (score = 19.23) (Supplementary Table 3) [36]. Univariate cox regression analysis indicated 29 genes with potentially prognostic impact at p < 0.01 (Supplementary Table 4).

The expanded dataset of 3,108 CDEGs links cell-cycle, ECM remodelling and metabolic alterations with increasing disease stage

To increase coverage of the BLCA alterations, we further extracted Concordantly Differentially Expressed Genes (defined in Methods; CDEGs, n = 3,108) from the comparisons of cancerous stages to normal adjacent urothelium and investigated the differentially activated pathways they represent, using ssGSEA. Similar to the analysis of monotonicity, we sought to identify pathways with monotonically up- or down-regulated activation scores across T0 and BLCA stages. Results indicated gradually stronger activations of several mitotic processes, positive regulation of the canonical Wnt pathway, mTORC1 signaling, expression of MYC targets, degradation of anaphase inhibitors, metabolism of nucleotides, mobility of formins, and the TNFR2/non-canonical NF-kB pathway (Figure 3, Supplementary Table 5). Conversely, consistently reduced expression versus controls was observed in genes involved in lipid and fatty acid catabolic processes, in the metabolism of heme, in pathways promoting muscle cell differentiation and, interestingly, in genes regulating the circadian clock (Figure 3). Regulon activity per sample was estimated and respective scores between cancerous stages and normal tissue were compared. GATA3 and GLI2 were the only regulons that significantly associated with increasing malignancy with their activity being downregulated (Figure 3).

• Download figure
• Open in new tab

Figure 3:

Excerpt of the pathways showing a monotonal increase or decline in their activation scores with higher stage.

To further investigate genome wide alterations and coexpression patterns associated with cancer progression, an integrated network-pathway analysis of the different disease stages was performed, using the 3,108 CDEGs. Out of 9,659,664 weighted gene pair interactions, the top 5,600 were used as the most representative to build the corresponding networks, applying Louvain clustering to identify topologically connected communities (modules). Retrieved gene clusters have high degree of coexpression, and thus changes in their composition across stages reflect gains/losses in pathway rewiring. We used the Gene Ontology – Biological Process library to identify molecular processes affected by changes in coexpression inside the top 5 largest in size (based on number of genes) communities (Figure 4, Supplementary Tables 6-11). The analysis revealed large differences in gene coexpression between T0 and cancerous samples with 4 out of the 5 largest in size communities associating clearly with specific biological processes.

• Download figure
• Open in new tab

Figure 4:

Biological process analysis of the largest in size coexpressed communities identified in each BLCA stage, using the 3,108 CDEGs between T0 and cancerous stages. A) Coherent communities identified and characterized across T0 and disease stages. Numbers in parentheses show the fraction of genes with existing Biological Process annotation with respect to the total number of genes found to be coexpressed in each of the community B) Barplots of the most significantly enriched biological processes per community depicting number of coexpressed genes for each. C) Hub genes identified across the studied conditions based on the betweenness centrality scores (y axis).

Three out of the top 5 communities were consistently detected in all BLCA stage networks. Based on examination of their enriched processes, these were labeled as 1) the cell-cycle community, 2) the ECM and developmental community 3) the metabolic and translational community (Figures 4A, 4B)

When investigating the cell cycle community and specifically the processes of G1/S transition of mitotic cell cycle, cell-cycle G2/M phase transition, regulation of cell cycle phase transition, regulation of chromosome organization and nuclear division, most of the participating genes were found to be coexpressed in all BLCA stages (Figure 4B). However, the gene size of the cell cycle community increased with higher stage (Ta n = 118, T1 n = 148 and MIBC n = 168-170), and interestingly, the proportion of the genes with unknown or unclear function (lacking GO annotation) also increased (11.9% for Ta, 29.8% for T4). Among the genes with cell cycle annotation (n=184), 80 genes were coexpressed in all cancerous stages and were also upregulated in malignancy compared to T0, possibly forming the backbone of cell cycle progression in BLCA (Supplementary Table 13). However, only 36 of them had a monotonal increase in the expression levels with higher stage (Supplementary Table 13). To detect potential drivers of variation in the cell cycle co-expression profiles of stages, we looked into the betweenness centrality scores of the genes inside each of the stage cell cycle communities (Figure 4C), and found CDCA5 and KIF2C to be hub genes in Ta tumors, FOXM1 and AURKB in T1, CDT1 and SMC4 in T2, CCNB1 and RRM2 in T3 and KIF14 in T4 tumors (Figure 4C).

The community of ECM and developmental processes was enriched in cell-cell communication and cell-matrix interactions, in responses to microenvironmental stress, as well as in differentiation programs of epithelial, mesenchymal and stem cells (Figure 5B). The biological process of extracellular matrix organization included coexpressions of 15-36 genes of which COL13A1, FGFR4, FOXF2 and SCUBE1 were coexpressed only in T0 compared to disease stages, whereas 26 genes (including COL6A1/A2, COL16A1, MFAP5, MMP11) were coexpressed in malignancy but not in T0. In line with recent observations [37], we noticed that T0 presented with an active ECM remodeling profile. Sixteen of the ECM associated genes were coexpressed both in T0 and in the NMIBC stages, including genes promoting ECM degradation and tumor cell migration (MMP2, CTSK, PDPN), fibrotic collagens (COL1A2, COL6A3, COL14A1, COL15A1), and pro-angiogenic factors (PDGFRA, RECK). Notably, coexpression in the T0 samples was predicted to be driven by the cancer stem cell marker ALDH1A2 and MFAP4 (Figure 4C), while in Ta tumors, by COL16A1 and CLIP3. Since CLIP3 interacts with both AKT1 and AKT2 [38], CLIP3 may have an important role in the early AKT/PI3K/mTOR axis of carcinogenesis. For the same process of extracellular matrix organization, compared to Ta, T1 tumors had gains in ADAMTSL2, while compared to T1, T2 tumors had gains in MMP11 and LRP1. Similarly compared to T2, T3 tumors had coexpression gains in ITGA10 and in genes of the endocytic pathway (ABL1, CYP1B1) whereas compared to T3, T4 tumors had gains in DDR2 and JAM2.

The community of metabolism and translation encompassed mitochondrial, translational and multiple metabolic processes being activated during carcinogenesis, and was more profound in the T1 and more advanced tumors. The processes of cellular respiration, translational initiation, mRNA catabolic process, nonsense mediated decay and protein targeting to ER were consistently enriched in most of BLCA stages. The results highlighted a set of 12 genes commonly co-expressed across stages for these processes, including COX7B, DLD, NDUFS4, UQCRFS1, PAIP2, RPL15, RPL30, RPL7, RPS23, RPS27, RPS27A, RPS4X, with RPL36A and YTHDF3 being T2-specific, COQ10B, EIF4E3, MTIF2, NCPB1 as T3-specific while EIF4E1, ISCU and GSPT2 were T4-specific.

Besides the abovementioned consistently detected communities in BLCA, a community enriched in processes of immune cell activation and cytokine secretion, was identified in the T0 and the MIBC stages and involved both innate and adaptive responses, as well as, processes of immune cell adhesion and migration. It was the least variable community with respect to different processes and respective genes coexpressed across the different stages. Strikingly, T0 and the MIBC stages had similar profiles. Out of the 17 genes of the process of T-cell activation that were commonly coexpressed at the MIBC stages, 15 were also coexpressed in the T0 samples. To further investigate these observations, the transcriptome data per sample were deconvoluted into relative abundances of immune cell populations using CIBERSORT, and cell fractions between disease stages were compared. Significant results were obtained for the following populations: CD8+, activated CD4+, activated NK, Monocytes, Macrophages M2 and activated Dendritic cells (Supplementary Figure 4). Results indicated differential commitment of immune cells to T0 and BLCA stages. T0 (n = 37) were significantly more infiltrated with CD8+ (p = 0.046) and monocytes (p = 4.6E-4) than tumor samples (n = 313), verifying the presence of the immune community in the T0 samples. However, compared to tumor, T0 samples had significantly less abundance of activated CD4+ cells (p = 5.28E-3), of macrophages (p = 16.8E-5), of activated dendritic cells (p = 0.0002) and of activated natural killer cells (p = 0.015). Generally, NMIBC had lower infiltration burden than MIBC, while compared to other BLCA stages, Ta tumors (n = 34) were significantly more infiltrated with activated dendritic cells (p = 0.024). Abundance of CD8+, of activated NK cells, and of M2 macrophages increased linearly with higher malignancy. Interestingly, AIF1 a gene that promotes macrophage survival [39], was found to be driver of immune coexpression in the T4 tumors.

RNA-seq independent analysis validates presence of communities and prognostic genes

In lack of another dataset comprising all the disease stage spectrum of primary BLCA incidents, the observed alterations in the discovery set were investigated for their reproducibility in the TCGA-BLCA publication [9]. Although the validation set is not perfectly suited since it includes only stages T2-T3-T4, we used it due to the RNA-seq technology which offers a good reliability at the validated findings. Differential expression analysis among the available stage comparisons (T3 vs T2 and T4 vs T3) in the TCGA data validated 74 of the 157 monotonal genes in the comparison T3vsT2, and none in the comparison T4vsT3 (Supplementary Table 13). Cox regression analysis of the 157 genes in the TCGA data validated 12 genes having a prognostic value (Supplementary Table 4), including IMPDH1, MRPL37, MED19, ENO1, ANLN, GTPBP4, MLF2, higher levels of which associating with worse survival and CBX7, ZFP2, AKAP7, ICA1L, CDC14B higher levels of which associating with better survival probability. To validate as possible the coexpression analysis findings, stage specific coexpression networks were also created using the TCGA data, and were clustered with the Louvain algorithm. GO-Biological process analysis of the communities validated the differential segregation of the cell-cycle, extracellular matrix and immune activation processes to distinct communities (Supplementary Figure 5).