Urothelial cancer gene regulatory networks inferred from large-scale RNAseq, Bead and Oligo gene expression data.

BACKGROUND: Urothelial pathogenesis is a complex process driven by an underlying network of interconnected genes. The identification of novel genomic target regions and gene targets that drive urothelial carcinogenesis is crucial in order to improve our current limited understanding of urothelial cancer (UC) on the molecular level. The inference of genome-wide gene regulatory networks (GRN) from large-scale gene expression data provides a promising approach for a detailed investigation of the underlying network structure associated to urothelial carcinogenesis.
METHODS: In our study we inferred and compared three GRNs by the application of the BC3Net inference algorithm to large-scale transitional cell carcinoma gene expression data sets from Illumina RNAseq (179 samples), Illumina Bead arrays (165 samples) and Affymetrix Oligo microarrays (188 samples). We investigated the structural and functional properties of GRNs for the identification of molecular targets associated to urothelial cancer.
RESULTS: We found that the urothelial cancer (UC) GRNs show a significant enrichment of subnetworks that are associated with known cancer hallmarks including cell cycle, immune response, signaling, differentiation and translation. Interestingly, the most prominent subnetworks of co-located genes were found on chromosome regions 5q31.3 (RNAseq), 8q24.3 (Oligo) and 1q23.3 (Bead), which all represent known genomic regions frequently deregulated or aberated in urothelial cancer and other cancer types. Furthermore, the identified hub genes of the individual GRNs, e.g., HID1/DMC1 (tumor development), RNF17/TDRD4 (cancer antigen) and CYP4A11 (angiogenesis/ metastasis) are known cancer associated markers. The GRNs were highly dataset specific on the interaction level between individual genes, but showed large similarities on the biological function level represented by subnetworks. Remarkably, the RNAseq UC GRN showed twice the proportion of significant functional subnetworks. Based on our analysis of inferential and experimental networks the Bead UC GRN showed the lowest performance compared to the RNAseq and Oligo UC GRNs.
CONCLUSION: To our knowledge, this is the first study investigating genome-scale UC GRNs. RNAseq based gene expression data is the data platform of choice for a GRN inference. Our study offers new avenues for the identification of novel putative diagnostic targets for subsequent studies in bladder tumors.

Background

Urothelial cancer (UC) is a heterogeneous disease with risk factors that include smoking, contact to chemicals and age [1]. Urothelial tumors originate from the epithelial lining of the bladder and can progress from non-invasive to more aggressive muscle-invasive subtypes which penetrate the deeper tissue layers of the bladder. The non-invasive tumor stages can be treated by transurethral resection, chemo- and intravesical therapy, whereas for invasive stages cystectomy, radiotherapy and chemotherapy are preferred [1,2]. Monitoring of UC is very expensive as its recurrence rate is high [3]. An understanding of the mechanistic interplay between individual genes and proteins that drive the development and progression of UC is therefore a high priority. System-based approaches allow us to investigate the underlying network structure associated with carcinogenesis and thus facilitate a novel perspective for the identification of molecular targets that drive urothelial carcinogenesis. The inference of gene regulatory networks (GRN) from large-scale gene expression data of tumor samples from various grades and stages is a promising approach for the identification of novel putative targets in cancer [4-6].

A GRN is a mathematical description of the dependencies within a gene expression dataset. Currently, a large arsenal of gene regulatory network inference methods have been developed [7-9]. The most popular methods are based on mutual information which is a dependency measure that can be estimated for all pairs of genes in a gene expression dataset. In this study we infer GRNs by the application of the BC3Net algorithm which is based on the aggregation of an ensemble of C3Net gene regulatory networks [4]. The C3Net algorithm selects a maximum of one gene neighbor for each gene on the basis of the strongest mutual dependency that is statistically significant. For a gene expression dataset with n genes we thus obtain a sparse network with at most n interactions. The BC3Net generates an ensemble of C3Net networks from bootstrap datasets, i.e., by sampling a dataset with replacement, that are subsequently aggregated to a weighted network. We have reported that BC3Net was shown to produce biological meaningful results [4-6,10,11]. The hub genes of GRNs that were inferred from large-scale cancer gene expression data were observed to provide promising putative novel target genes for cancer such as G-protein coupled receptors and transmembrane proteins [10,12].

Methods

Preprocessing and sample information for the Illumina RNAseq, Bead array and Affymetrix Oligo microarray gene expression dataset

We preprocessed three large-scale urothelial cancer gene expression datasets from a) Illumina RNAseq (179 samples) [13], b) Illumina Bead array (165 samples) GSE13507 [14,15] and c) Affymetrix oligo microarray (188 samples) platform [16-19]. An overview of the tumor stage information for the three datasets is shown in Table 1. In this table we distinguish 6 tumor stages, namely, pTcis, pTa, pT1, pT2, pT3 and pT4. For each of these stages we list the number of available samples provided by the three platforms.

Table 1

Tumor stages across the three datasets

Tumor stage	RNAseq	Bead	Oligo
	(# samples)	(# samples)	(# samples)
pTcis			5
pTa	1	24	35
pT1	1	80	16
pT2	57	31	22
pT3	91	19	51
pT4	29	11	29
pT2-4			30
Total	179	165	188
NMI	2	104	56
MI	177	61	132
Genes	20,161	18,956	12,495

The RNAseq dataset consists mainly of muscle-invasive UC tumor samples. 30 samples of the Oligo dataset corresponded to the muscle invasive stages pT2 to pT4 and were not assigned a specific stage.

RNAseq gene expression dataset from TCGA

The RNAseq gene expression dataset was retrieved from The Cancer Genome Atlas bladder cancer TCGA project [http://cancergenome.nih.gov/] [13]. We used the preprocessed RNAseqV2 normalized count expression values based on RSEM (RNA-Seq by Expectation-Maximization) [20,21] as provided by TCGA and clinical information such as the TCGA barcode identifier, sample type and tumor histology by the bcr_aliquot_uuid identifiers. We extracted gene expression data of primary solid tumors for a total of 179 samples with histology stage information (march 2014). A total of 177 of the 179 selected tumor samples represent muscle invasive carcinoma stage pT2 or above. We performed a log-transformation log(1+p) on the count expression values. The resulting gene expression matrix consisted of 20,161 entrez genes and 179 samples. Genes with a zero standard deviation were removed from the dataset.

Illumina Bead array gene expression dataset

We used the processed matrix series Illumina Bead gene expression data from GEO GSE13507 [14,15]. The dataset comprises 257 samples from the Illumina human-6 v2.0 expression Beadchip microarray platform. The dataset consists of tumor samples from 62 muscle invasive, 104 non-muscle invasive, 23 recurrent non-invasive bladder cancer and samples that we excluded for the network inference representing 58 mucosae surrounding cancer and 10 normal mucosae. We assigned Illumina identifiers to entrez gene id and gene symbols using the illuminaHumanv2.db annotation bioconductor package. From the 43,148 Illumina identifiers for a total of 20,481 an entrez identifier was available. The remaining 22,667 features were not considered for the analysis. In total, we selected 165 primary bladder cancer samples (Table 1).

Affymetrix Oligo microarray gene expression dataset

We used a third UC dataset from Affymetrix gene expression data comprising 183 samples from 4 different datasets. We extracted 93 (U133plus2) samples from GSE31684 [16], 46 (U133A) samples from GSE3167 [17], 30 (U133A) samples from GSE5287 [18] and 19 (U133A) samples from GSE37317 [19]. We considered only probe sets that were present in both array types U133a and U133plus2. We combined the U133plus2 samples and the U133a samples using the matchprobes bioconductor package [22]. We normalized the microarray samples using RMA and quantile normalization [23] using log2 expression intensities for each probe set. As a summary statistic for multiple probesets that match to the same entrez gene identifiers we used the median expression value. Entrez gene ID to Affymetrix probe set annotation was obtained from the hgu133plus2.db and hgu133a.db R package. We excluded all probe sets from our analysis that remained unmapped to entrez identifiers. The resulting expression dataset consisted of 12,495 genes and 188 samples.

BC3Net gene regulatory network inference

We inferred our bladder cancer GRN using C3Net and the “B”agging version of the C3Net [24,25] algorithm called BC3Net [4]. The BC3Net infers an ensemble of C3Net gene regulatory networks from bootstrap generated datasets that are subsequently aggregated to a weighted GRN. We defined an ensemble of B=100 independent bootstrap datasets that were generated from a given gene expression dataset D. For each bootstrap data set a GRN was inferred using C3Net [24,25]. Edges with non-significant mutual information values were subsequently rejected using a non-parametric test with a Bonferroni multiple testing correction for a significance level α=0.05. The null distribution of mutual information is generated from sample-gene label permutations of the original gene expression matrix. For the network inference we used a Pearson Estimator [8,26]

where ρ denotes the Pearson correlation coefficient. The inferred ensemble of GRNs was aggregated into a weighted network . The weights of the inferred interactions give the frequency how often an interaction was observed in the C3Net network ensemble and are denoted as ensemble consensus rate (ECR). For each inferred weighted edge in the network the statistical significance was estimated by a Binomial test. For multiple testing correction Bonferroni was used with a significance level α=0.05.

Relevance networks

For the inference of relevance networks [27] we used the WGCNA R Package [28] and the CLR [29] implementation provided in the minet R-Package [30]. Interactions were defined for WGCNA by hard thresholds on the absolute Pearson correlation matrix and for CLR by hard thresholds on the z-score transformed mutual information matrix that was estimated using a Pearson Estimator [8,26].

Cancer census genes

The Cancer Gene Census (CGC) [31] (version download 10-01-2014) [http://www.sanger.ac.uk/genetics/CGP/Census/] provides information about genes with somatic mutations that are associated to different types of cancer. We used the entrez identifiers of the defined cancer census genes.

Gene ontology gene sets

For our analysis, we obtained the Gene Ontology [32] annotation for entrez gene IDs from Bioconductor [22] annotation packages org.Hs.eg.db and GO.db.

Gene family gene sets

We retrieved gene family protein tag information and entrez identifiers of the genes in the HGNC database [http://www.genenames.org/genefamilies]. We defined gene family gene sets for groups of genes that shared the same HGNC protein family tag. From the HGNC database we gathered a total of 587 gene family gene sets comprising a total of 16,722 entrez genes.

Gene sets of co-localized adjacent genes

For the identification of genomic regions with enriched subnetworks of interacting genes we defined gene sets of genes that were adjacently located within a chromosomal region (co-located) from overlaping sliding windows along the human chromosomes. We defined gene sets from 1Mb (mega bases) sliding windows along the human chromosomes with a 500Kb (kilo bases) overlap between adjacent windows. The gene sets of co-located genes were defined for chromosome regions of 1Mb with 500Kb overlap to mimic the extend for co-expressed gene clusters [33].

Gene pair enrichment analysis (GPEA)

The GPEA facilitates the identification and ranking of significant subnetworks of defined gene sets for a given network. For p genes there is a total of N=p(p−1)/2 different gene pairs. If there are p genes for a particular gene set (S) then the total number of gene pairs for this gene set is m=p(p−1)/2. When a network G contains n interactions, of which k interactions are among genes from the given gene set S, then a p-value for the enrichment of gene pairs of this gene set S can be calculated from the following hypergeometric distribution

This p-value gives an estimate for the probability to observe k or more interactions between genes from a given gene set S.

We performed a GPEA analysis for the inferred GRNs for ∼ 8,000 gene sets of GO biological process (≥3 and <1000 genes), ∼ 500 gene sets of gene families (≥3 genes) and ∼ 4,000 gene sets of co-located genes (≥3 genes). For the analysis the inferred networks are expected to show a strong association to gene sets of a biological functional and spatial context. Therefore, we considered a more stringent significance level of α=0.001 (10−3) relative to the number of performed test in the range of 103. Further, we considered a Bonferroni multiple testing correction.

Network centrality measures

For the network analysis we measured the degree centrality and edge density [34]. The degree centrality was defined as the total number of direct neighbors of a gene g of an undirected gene regulatory network. The edge density of a network was the number of edges divided by the maximal number of possible edges. For an undirected network this number was given by n(n−1)/2, whereas n is the total number of genes.

Protein interaction databases

We gathered and processed interactions from Biogrid [35] (15,337 genes, 135,732 interactions; version biogrid.3.2.11), Intact [36] (10,029 genes, 63,968 interactions; version intact.230314), Mint [37] (7,106 genes, 26,834 interactions; version mint.2013-03-26), Hprd [38] (9,672 genes, 39,233; version hprd.072010), String [39] (20,770 genes, 4,850,628 interactions, version 9.1). Further, we considered the largest manually curated human signaling network [40] (6,306 genes and 57,090 interactions, version 6) which we denote in the text as SingNet (http://www.cancer-systemsbiology.org/dataandsoftware.htm), a pathway protein interaction network extracted from the bioconductor package graphite [41] (6,243 genes, 78,201 interactions; KEGG, Reactome, NCI and Spike) and the integrative network from ConsensusPathDB (CPDB) [42] (16,619 genes, 485,277 interactions; version Dec 2014). We assigned their entrez gene identifiers mapping when available from the interaction database or converted the identifiers (e.g. uniprot identifiers) to entrez identifiers using the annotation from the bioconductor package org.Hs.eg.db and uniprot database [43].

Quantitative comparison of experimental interactions in RNAseq, Bead and Oligo UC GRN

We used the interactions from the Biogrid, Intact, Mint, Hprd, CPDB, SingNet, graphite and String database seperately as global reference networks for the GRN and measured the number of true, false positive (TP, FP), true, false negatives (TN, FN) and F-score to compare the performance of the three inferred gene regulatory networks. The F-score measure gives a weighted average of the precision and recall .

For the local subnetwork based network inference performance comparisons we used the String network as a reference network. We compared the cumulative log transformed F-score distribution separately for commonly significant GO Biological Process, genomic co-located genes and gene family subnetworks between the Oligo, Bead and RNAseq UC GRNs.

For all subnetworks and pairwise network comparisons we performed a hypergeometric test for the number of shared interactions between two networks is not larger than expected by random chance. For the subnetwork analysis we consider FDR multiple testing correction [44].

Results

Urothelial cancer (UC) gene regulatory networks (GRN)

For the identification of molecular targets for UC from a network-based perspective we inferred BC3net GRNs from RNAseq, Bead and Oligo UC gene expression datasets. The giant connected component of the inferred RNAseq UC GRN consisted of 18,952 genes, the Bead UC GRN of 20,140 genes and the Oligo UC consisted of 12,492 genes (Figure 1A). In the following we compared the global network and local structural properties between the three networks.

Figure 1

Overview of the comparisons between the inferred RNAseq (A), Bead (C) and Oligo (E) UC GRN. The pairwise comparisons of the three networks are shown in (B, D, F) and triple comparison in (G). We compared the networks on the gene interaction level and gene set subnetwork level for Gene Ontology biological processes, gene families and co-located gene sets of 1 Mb genomic regions.

The global network properties among the three networks were highly similar. The edge density (d∼0.001) of the Oligo has a slightly higher edge density compared to the RNAseq and Bead UC GRN. The degree distribution of the UC GRN follow a power law distribution with exponents αRNAseq=4.09, αBead=4.16αOligo=3.73. The average shortest path length for all three networks was pRNAseq=4.17, pBead=4.45, pOligo=4.43 genes (measured with the Dijkstra distance [45]).

On the local structural level of individual interactions we observed that the three networks are highly dissimilar. We performed a pairwise comparison and joint comparison of all three networks to quantify the number of shared interactions. The percentages of shared interactions were quantified from the union of all interactions of two networks and for the joint comparison from the union of all interactions of the three networks. The GRN networks shared only a total of 6,676 (2.17%, RNAseq/Bead), 5,321 (2.04%, RNAseq/Oligo) and 2,648 (1.22%, Bead/Oligo) interactions which corresponded to subnetworks among 3,312 to 7,361 genes. In total, we found that only 1,002 (0.25%) interactions were shared across the three GRNs and corresponded to a subnetwork among 1,539 genes. An overview of our gene expression data and inferred gene regulatory networks on the interaction level is shown in Figure 1.

Functional analysis of the inferred Oligo, Bead and RNAseq UC GRN

In this section we highlighted the key biological processes of the three UC GRNs and their association to known cancer genes and performed a comparative analysis between inferred GRN, relevance and PPI networks. We identified the most prominent subnetworks for known biological processes of the inferred Oligo, Bead and RNAseq UC GRN by a functional enrichment analysis for gene pairs (GPEA). The association of the identified biological processes to known cancer hallmarks was quantified by a subsequent enrichment analysis of cancer census genes.

The GPEA analysis was performed for Gene Onotology (GO) biological process for all terms with ≥3 and <1000 genes. For the RNAseq UC GRN we observed a total of 10.3% significant GO terms. In contrast, for the Bead and Oligo UC GRN only 4.94% and 5.38% of all tested GO terms were significant Figure 1. For all GRN networks we observed that 50% of the identified significant GO terms were also enriched by cancer census genes. A total of 91% (RNAseq), 88% (Bead) and 93% (Oligo) of the cancer genes were present in the selected set of significant Gene Ontology biological processes.

From all significant GO terms that we identified 299 (∼55.4%) GO terms were common across the three UC GRNs. We observed a wide variety of common biological processes with a pronounced representation of immune related processes, cell cycle, catabolic processes such as proteolysis, chromatin organization, metabolism, cell adhesion, cell migration, cell differentiation and development including keratinization and angiogenesis. A complete list of the significant terms for the individual analyses is given in the Additional file 1: Tables S1, S2 and S3. An overview of the functional landscape of the common significant terms among the GRN networks is shown in Additional file 1: Figure S1.

In order to evaluate our results we compared the fraction of cancer associated biological processes between the BC3Net and C3Net GRN, WGCNA and CLR relevance networks and PPI networks from graphite, SingNet and CPDB. The analysis was performed separately for the Oligo, Bead and RNAseq gene expression data. For the analysis we generated relevance networks by hard thresholds for 0.1 to 0.9 percentiles of the absolute correlation matrix (WGCNA) and the z-score transformed matrix (CLR). C3Net inferred interactions were weighted by the respective mutual information value and for BC3Net by the ensemble consensus rate (ECR). For C3Net and BC3Net the analysis was performed on the entire network and for an ensemble of hard thresholds ranging from 0.1 to 0.9 percentiles. Figure 2 shows the fraction of biological processes that were identified from the GPEA analysis (α=0.001, Bonferroni) with a significant enrichment of cancer census genes. For all 3 datasets the BC3Net showed the largest fraction of cancer associated significant biological processes (∼50%). CLR and WGCNA showed a low performance on the Oligo gene expression dataset (25% to 35%) that is comparable to PPI networks. We also observed that CLR shows a prominently improved performance compared to WGCNA.

Figure 2

GPEA analysis of GO biological process gene sets for relevance networks, C3Net and BC3Net GRNs. Shown are the fraction of significant GO biological process terms with significant enrichment of cancer census genes.

Gene pair enrichment analysis of gene subnetworks of co-located genes

Gene expression profiles that are influenced by genomic and epigenomic alterations can elucidate dependency structures of co-located genes and link to novel genomic target regions which are specific to urothelial cancer. In order to identify genomic cancer target regions with significant subnetworks in the RNAseq, Bead and Oligo GRNs we perfromed an enrichment analysis for gene pairs in gene sets from genome-wide 1 Mb genomic regions of co-located genes.

We observed 11.68% significantly co-located gene subnetworks for the RNAseq UC GRN. In contrast, for the Bead and Oligo UC GRN we identified only 2−3%. Figure 3 shows the most prominent co-located gene subnetwork for the GRNs of the RNAseq, Bead array and Oligo UC dataset. For the three GRNs the top 50 chromosomal regions with a significant GRN subnetwork are shown in Tables 2, 3 and 4 (for full tables see Additional file 1: Tables S4, S5 and S6). We reviewed the literature for the most prominent identified genomic region and their association to UC for each GRN. For the RNA-seq UC GRN the most prominent gene subnetwork was located on chromosome locus 5q31.3 and represents a protocadherin gene cluster. In UC, the loci 5q31.2−q32 has been associated with losses in a low fraction of UC tumors [46,47]. In [48,49] an epigenetic analysis was performed on free DNA derived from blood serum samples from UC patients of the protocadherin PCDH10 an PCDH8. The studies of [48,49] showed that the methylation patterns of PCDH10 and PCDH8 were significantly associated with stage, grade, recurrence and tumor size. The 5q31.3 locus was also described in Wilms tumor to be epigentically silenced [50]. For the Bead UC GRN the most prominent co-located gene subnetwork was located on chromosome locus 8q24.3. The 8q24.3 locus is a common gain loci in multiple cancers and was also confirmed in multiple UC cell lines [51]. For the Oligo UC GRN the chromosome locus 1q23.3 contained the most prominent co-located gene subnetwork. In [52] a gain of 1q23.3 was identified from free DNA in urine samples from UC patients. Further, [52] validated a selected candidate PFND2 that is located in 1q23.3 in an independent set of urothelial cancer tumors. PFND2 was significantly amplified and overexpressed and showed association to increasing stage and tumor grade.

Figure 3

Significant genomic UC GRN subnetworks defined from co-located genes of 1 Mb genomic regions.(A) Overview of the GPEA analysis showing all significant genomic subnetworks for the RNAseq (yellow), Bead (blue) and Oligo (red) UC GRNs. (B-D) Shown are the most prominent genomic subnetworks derived from the RNAseq (B), Bead (C) and Oligo (D) UC GRNs.

Table 2

GPEA analysis of 1 Mb genomic regions gene sets for the RNAseq UC GRN

chr	Locus	Start	Size	Edges	p-value	Census
chr5	q31.3	140000001	74	159	3.6673e-222
chr17	q21.2	39000001	61	136	7.6570e-204
chr17	q21.2	38500001	49	119	1.1284e-194	RARA, SMARCE1
chr6	p22.2	25500001	50	105	1.0881e-163
chr6	p22.2/p22.1	26000001	48	98	2.9016e-153
chr5	q31.3	140500001	54	103	1.6926e-152
chr19	q13.43	57500001	41	90	3.9091e-150
chr8	q24.3	145000001	43	81	1.5314e-127	RECQL4
chr19	q13.43	58000001	46	75	8.4630e-111
chr16	p11.2	30000001	54	80	4.0300e-109
chr21	q22.11	31500001	38	67	5.2390e-107
chr1	q21.3	152500001	43	69	2.0553e-103
chr9	q34.3	139500001	67	83	7.6570e-99
chr17	q25.3	79500001	47	67	3.9091e-94	ASPSCR1
chr19	q13.31	44000001	37	60	4.0300e-94
chr16	p11.2	30500001	49	68	1.7329e-93	FUS
chr8	q24.3	145500001	37	58	4.8360e-90	RECQL4
chr11	q13.1/q13.2	65000001	46	64	8.8660e-90
chr3	p25.3	9500001	34	55	3.1837e-88	FANCD2, VHL
chr5	q31.3	139500001	42	60	3.1434e-87
chr6	p22.1	27000001	32	51	8.8660e-83	HIST1H4I
chr1	q21.3	152000001	34	52	5.2390e-82
chr16	p11.2	29500001	45	59	1.1687e-81
chrX	p11.23	48500001	50	62	1.8941e-81	GATA1, TFE3, WAS
chr6	p21.33/p21.32	31500001	79	78	1.2896e-79
chr19	p13.3	500001	51	60	7.6570e-77	FSTL3, STK11
chr1	p36.33	1000001	48	58	1.7732e-76
chr19	p13.3	1000001	43	54	1.1687e-74	STK11, TCF3
chr1	q21.3	150500001	35	48	1.5314e-72	MLLT11, ARNT
chr1	q22	155000001	42	51	3.3852e-70	MUC1
chr8	q24.3	144500001	39	49	9.2690e-70
chr19	q13.41/q13.42	53000001	31	44	1.6926e-69
chr19	p13.11/p12	19000001	33	45	5.2390e-69
chr16	p13.3	1	52	56	5.6420e-69	AXIN1
chrX	p11.23/p11.22	49000001	44	51	4.4330e-68
chr11	q13.1/q13.2	65500001	46	52	7.6570e-68
chr6	p22.1	27500001	36	46	2.1359e-67
chr19	q13.33	49500001	70	64	2.9822e-66
chr1	q23.3	161000001	36	45	1.8941e-65	FCGR2B, SDHC
chr19	q13.31/q13.32	44500001	36	45	1.8941e-65	CBLC, BCL3
chr3	p21.31	48500001	42	48	1.1284e-64	NCKIPSD
chr19	q13.43	58500001	30	41	1.4911e-64
chr16	p13.3	1500001	59	57	2.1762e-64	TSC2, TRAF7
chr1	q21.3	151000001	37	45	2.4583e-64	MLLT11
chr11	p15.5	1	50	52	5.2390e-64	HRAS
chr1	q21.3/q22	154500001	42	47	7.6570e-63	MUC1
chr19	q13.12	36000001	48	50	1.9747e-62
chr11	q13.1	64500001	56	54	4.4330e-62	MEN1
chr19	p13.2	11500001	30	39	1.9747e-60
chr11	q13.2	67000001	36	42	1.1284e-59

For each significant genomic region the chromosome (chr), chromosome band (locus), nucleotide base start site of the genomic region (start), number of genes of the gene set (size), number of edges of the significant subnetwork (edges), Bonferroni adjusted p-value of the subnetwork (p-value) and a list of genes in the significant subnetwork that are present in the cancer census (census).

Table 3

GPEA analysis of 1 Mb genomic regions gene sets for the Bead UC GRN

chr	Locus	Start	Size	Edges	p-value	Census
chr8	q24.3	145000001	43	60	2.08131e-90	RECQL4
chr6	p22.2/p22.1	26000001	48	54	2.23839e-73
chr6	p22.2	25500001	50	54	2.00277e-71
chr1	q21.3	152500001	43	46	1.76715e-63
chr17	q25.3	79500001	47	41	4.31970e-51	ASPSCR1
chr6	p22.1	27000001	32	34	7.06860e-51	HIST1H4I
chr8	q24.3	145500001	37	35	2.51328e-48	RECQL4
chr8	q24.3	144500001	39	30	1.02102e-37
chr8	p11.23/p11.22	37500001	20	21	3.49503e-35	WHSC1L1
chr1	q23.3	161000001	36	27	2.04204e-34	FCGR2B, SDHC
chr11	q13.2	67000001	36	25	7.06860e-31
chr17	q21.32	46000001	31	23	1.80642e-30
chr17	q21.32/q21.33	46500001	31	23	1.80642e-30
chr1	q21.3	153000001	43	27	3.29868e-30
chr11	p15.5	1	50	29	1.06029e-29	HRAS
chr16	p13.3	1	52	29	1.02102e-28	AXIN1
chr19	q13.43	57500001	41	25	5.10510e-28
chr1	q23.3	160500001	34	22	7.85400e-27	SDHC
chr1	p34.3/p34.2	40000001	23	18	2.74890e-26	MYCL1
chr6	p22.1	27500001	36	22	1.02102e-25
chr1	p34.2	40500001	22	17	5.89050e-25
chr1	q21.3	150500001	35	20	7.46130e-23	MLLT11, ARNT
chr1	q21.3/q22	154500001	42	22	9.42480e-23	MUC1
chr6	p21.32/p21.31	32500001	42	22	9.42480e-23	DAXX
chr9	q34.3	139500001	67	28	2.98452e-21
chr19	q13.43	58000001	46	22	4.71240e-21
chr1	q21.3	151000001	37	19	2.82744e-20	MLLT11
chr16	p13.3	500001	48	22	3.10233e-20
chr3	p21.31	49500001	38	19	7.85400e-20
chr1	q22	155000001	42	20	1.09956e-19	MUC1
chr11	q13.2	66500001	35	18	1.49226e-19
chr5	q31.3	140000001	74	28	5.89050e-19
chr22	q13.33	50000001	30	16	1.68861e-18
chr4	q13.2/q13.3	69500001	9	9	2.67036e-17
chr6	p21.1	42500001	33	16	3.69138e-17
chr12	q15	69000001	15	11	3.76992e-17	MDM2
chr12	q13.3/q14.1	57500001	37	17	3.92700e-17	CDK4, DDIT3
chr17	q11.2	26500001	49	20	4.71240e-17
chr19	p13.2	7500001	42	18	1.02102e-16
chr9	q34.3	140000001	35	16	2.43474e-16
chr20	q13.12	43500001	39	17	2.43474e-16	SDC4
chr22	q13.33	50500001	32	15	5.89050e-16
chr5	q31.3	140500001	54	20	2.04204e-15
chr1	q21.3	152000001	34	15	3.76992e-15
chr7	p15.2/p15.1	27000001	22	12	3.76992e-15	JAZF1, HOXA11, HOXA13
chr16	q12.2/q13/q21	56500001	34	15	3.76992e-15	HERPUD1
chr2	q35	219500001	47	18	5.49780e-15	FEV
chr17	p13.1	7000001	68	23	5.89050e-15	TP53
chr11	q13.1	64500001	56	20	8.24670e-15	MEN1
chr8	p11.21	41500001	12	9	1.02102e-14	KAT6A

For each significant genomic region the chromosome (chr), chromosome band (locus), nucleotide base start site of the genomic region (start), number of genes of the gene set (size), number of edges of the significant subnetwork (edges), Bonferroni adjusted p-value of the subnetwork (p-value) and a list of genes in the significant subnetwork that are present in the cancer census (census).

Table 4

GPEA analysis of 1 Mb genomic regions gene sets for the Oligo urothelial cancer GRN

chr	Locus	Start	Size	Edges	p-value	Census
chr1	q23.3	161000001	36	27	4.82850e-30	FCGR2B, SDHC
chr17	q25.3	79500001	47	30	5.15040e-28	ASPSCR1
chr1	q21.3	150500001	35	25	1.80264e-27	MLLT11, ARNT
chr11	p15.5	1	50	25	8.36940e-20	HRAS
chr8	q24.3	145000001	43	22	8.04750e-19	RECQL4
chr1	p34.3/p34.2	40000001	23	15	6.11610e-18	MYCL1
chr1	q23.3	160500001	34	18	3.86280e-17	SDHC
chr19	q13.2/q13.31	43000001	20	13	4.18470e-16
chr9	p21.1/p13.3	32500001	17	12	4.50660e-16
chr16	p13.2	8000001	7	8	1.41636e-15
chr6	p21.1	42500001	33	16	1.25541e-14
chr1	q21.2/q21.3	150000001	27	14	3.54090e-14	ARNT
chr8	p11.21	42000001	17	11	4.18470e-14
chr16	p13.2	8500001	9	8	2.06016e-13
chr16	q21/q22.1	66500001	42	17	1.31979e-12	CBFB
chr17	q12/q21.1/q21.2	37500001	31	14	1.67388e-12	ERBB2, CDK12, RARA
chr4	q13.3	74500001	14	9	5.47230e-12
chr6	p21.1	43000001	26	12	1.80264e-11
chr1	q21.3/q22	154500001	42	16	2.44644e-11	MUC1
chr4	q13.3	74000001	16	9	7.40370e-11
chr8	p11.23/p11.22	37500001	20	10	9.65700e-11	FGFR1, WHSC1L1
chr11	q13.2	67000001	36	14	1.03008e-10
chr9	p21.3	20500001	24	11	1.06227e-10	MLLT3
chr1	p34.2	40500001	22	10	6.75990e-10
chr11	p11.2	47000001	22	10	6.75990e-10	DDB2
chr13	q34	113500001	18	9	6.75990e-10
chr2	p22.1	39000001	11	7	1.22322e-09
chr13	q14.2	48500001	11	7	1.22322e-09	RB1
chr19	p13.12	15000001	27	11	1.41636e-09	BRD4
chr17	q25.3	80000001	23	10	1.67388e-09
chr17	p13.3	500001	19	9	1.83483e-09	YWHAE
chr1	p34.3	37500001	20	9	4.82850e-09
chr1	q22	155000001	42	14	6.75990e-09	MUC1
chr9	p21.3	21000001	25	10	9.01320e-09	CDKN2A
chr3	p25.3	9500001	34	12	1.06227e-08	VHL
chr16	q12.2/q13/q21	56500001	34	12	1.06227e-08	HERPUD1
chr19	q13.33/q13.41	51000001	49	15	3.21900e-08	KLK2
chr1	p34.3	38000001	18	8	3.86280e-08
chr11	q13.1/q13.2	65000001	46	14	7.40370e-08
chr11	q13.1/q13.2	65500001	46	14	7.40370e-08
chr19	q13.31	44000001	37	12	7.72560e-08
chr20	p13	3000001	28	10	8.69130e-08
chr6	p21.32/p21.31	32500001	42	13	9.97890e-08	DAXX
chr12	q15	69000001	15	7	1.31979e-07	MDM2
chr20	q11.22/q11.23	33500001	24	9	1.35198e-07
chr1	q21.3	153000001	43	13	1.80264e-07
chr9	p21.1/p13.3	33000001	20	8	2.12454e-07
chr22	q11.21	20500001	20	8	2.12454e-07
chr8	q24.3	144500001	39	12	2.60739e-07
chr11	q13.2	66500001	35	11	3.86280e-07

For each significant genomic region the chromosome (chr), chromosome band (locus), nucleotide base start site of the genomic region (start), number of genes of the gene set (size), number of edges of the significant subnetwork (edges), Bonferroni adjusted p-value of the subnetwork (p-value) and a list of genes in the significant subnetwork that are present in the cancer census (census).

Gene pair enrichment analysis of gene family subnetworks

A gene family is a group of duplicated genes with similar biological functions or biochemical activities which often form gene clusters of genes with chromosomal co-location. In this section we performed a GPEA for gene family gene sets and compared the results between the UC RNAseq, Bead and Oligo GRNs.

We found a total of 20.44% (103) significant gene family subnetworks in the RNAseq GRN (Additional file 1: Table S7), but only 8.2% (40) for the Bead UC GRN (Additional file 1: Table S8) and 9.1% (39) for the Oligo UC GRN were significantly enriched (Additional file 1: Table S9). However, we found a high agreement of the gene family subnetworks between our GRNs with a total pairwise overlap of 30−60% and among the three networks 25% of all identified gene family subnetworks (Figure 1D). It is noteworthy that 97.5% of the gene families identified by the Bead and 82% of the gene families identified by the Oligo network were also significantly enriched in the RNAseq GRN. In total 28 gene families were identified across all three UC GRNs which described CD molecules, keratin proteins (KRT), protocadherins (PCDHC), kalikrein proteins (KLK), zinc-finger transcription factors, metallothioneins and the immunoglobulin superfamily. An overview of the significant gene families that were common across all three GRNs are shown in Table 5.

Table 5

GPEA analysis of gene family gene sets for the UC RNAseq, Bead and Oligo GRN

Tag	Name	RNAseq	Bead	Oligo
CD	CD molecules	380/591/0	365/319/1.1e-140	351/319/8.9e-106
ZKRAB	Zinc fingers, C2H2-type with KRAB domain	338/764/0	306/289/5.3e-155	143/108/5.6e-64
RPL	L ribosomal proteins	59/119/1.2e-175	55/49/2.2e-59	44/72/3.4e-102
ZNF	Zinc fingers, C2H2-type	697/1307/0	641/515/4e-114	362/242/2.3e-53
HIST	Histones/Replication-dependent	67/236/0	61/154/5.7e-252	26/32/1.3e-48
HLA	Histocompatibility complex	24/46/1.1e-85	25/33/2.9e-57	19/22/9.1e-36
C1SET	Immunoglobulin superfamily/C1-set domain containing	38/51/1e-75	38/34/1.5e-46	37/26/1.1e-28
KLK	Kallikreins	17/26/1.2e-49	16/15/2.8e-26	14/14/1.7e-23
MT	Metallothioneins	14/12/3.5e-20	12/12/7.1e-23	10/10/3.5e-18
PCDHC	Cadherins/Protocadherins : Clustered	57/150/7.8e-242	56/24/6.6e-21	22/14/1.3e-17
IGD	Immunoglobulin superfamily/Immunoglobulin-like domain containing	233/175/3.6e-85	221/77/5.1e-22	177/67/2.6e-17
KRT	Keratins	55/87/1.4e-121	51/20/2.8e-17	35/20/1.5e-20
HOXL	Homeoboxes/ANTP class : HOXL subclass	52/54/1.6e-66	50/39/1.8e-46	43/20/5.2e-17
COLLAGEN	Collagens	46/37/3.4e-43	43/21/1.1e-21	33/17/5.5e-17
RPS	S ribosomal proteins	34/40/7.2e-59	32/14/3.1e-15	29/28/6.1e-38
PSM	Proteasome (prosome, macropain) subunits	45/25/4.1e-25	43/15/5e-13	42/26/7.7e-26
RBM	RNA binding motif (RRM) containing	209/114/2.6e-46	186/50/1.6e-12	151/72/2.3e-28
ENDOLIG	Endogenous ligands	230/77/7.9e-16	221/62/3.3e-13	192/62/2.8e-11
VSET	Immunoglobulin superfamily/V-set domain containing	161/95/1.4e-50	150/42/5.2e-14	110/33/9e-11
comI	Mitochondrial respiratory chain complex/Complex I	38/14/3.9e-12	38/12/2.6e-10	31/15/7.3e-15
UGT	UDP glucuronosyltransferases	20/15/4.8e-22	20/17/2.3e-27	7/5/1.4e-08
IFF2	Intermediate filaments type II, keratins (basic)	26/23/2.6e-33	24/8/3e-08	15/7/1.8e-08
comIV	Mitochondrial respiratory chain complex/Complex IV	16/13/8.5e-21	15/6/1.3e-07	12/7/6.3e-10
S100	S100 calcium binding proteins	21/12/1.1e-15	21/7/1.8e-07	17/7/1.1e-07
LNCRNA	Long non-coding RNAs	548/407/1.2e-78	468/168/2e-14	107/25/5.4e-06
complement	Complement system	35/18/3.8e-19	33/7/1e-04	30/9/9.6e-07
CYP	Cytochrome P450s	62/37/1.1e-33	56/10/0.00019	48/22/1.2e-17
SERPIN	Serine (or cysteine) peptidase inhibitors	36/21/2.1e-23	35/8/1.2e-05	30/7/0.00031

Shown are the commonly significant gene family subnetworks of the RNAseq, Bead and Oligo UC GRN, number of genes of the gene family subnetwork, the number of interactions and the Bonferroni adjusted P-value (Genes/Interactions/P-value).

Bladder cancer GRN degree centrality and hub genes

The identification of highly interactive central genes, i.e., hub genes of inferred and experimental network can provide promising targets for urothelial cancer. In this section we described individual hub genes of the gene regulatory network and review their functional role and relevance for the study of UC.

In order to compare the global structural properties of individual genes we performed a pairwise comparison of the degree centrality for 11,700 genes that are present among the three networks. The pairwise comparisons of the gene degree centrality across the three networks showed only a weak correlation. The degree ranks showed a slightly higher correlation between the RNAseq and Bead GRN of r=0.22 (RNAseq-Bead, p≤2.2e−16) compared to r=0.16 (Bead-Oligo, p≤2.2e−16) and r=0.16 (RNAseq-Oligo, p≤2.2e−16).

Hub genes of gene regulatory networks were observed to be highly dataset specific. Table 6 A, B, C shows the six most frequently observed hub genes for each of the inferred UC GRNs. In the following we describe the hub genes for which there is strong evidence for their relevance to cancer related properties. For example, the transmembrane protein HID1 that was observed as a major hubgene in the RNAseq GRN is reported to be downregulated in multiple cancers [53]; FER1L4 is a lncRNA reported to be prominently downregulated in gastric cancer [54], TTLL3 is described as a candidate cancer gene [55], RIF1 has been described to have anti-apoptotic properties in DNA repair [56] and SBNO1 (strawberry notch homolog 1) was reported in lung cancer [57]. For the UC Bead GRN, RNF17 (TDRD4) is a potential liver cancer CT antigen [58] and TMED7 was observed to be upregulated in a nasopharyngeal carcinoma cell line and described to act as a major immune system switch [59]. The Oligo GRN hubgene CYP4A11 was shown to promote angiogenesis and metastasis in lung cancer [60] and SLC38A3 (SNAT3) is a glutamine transporter and has been described as a marker for malignant glioma [61].

Table 6

The six major hub genes for the RNAseq, Bead and Oligo UC GRN

A) RNAseq UC GRN
Gene	Degree	Locus	Description/GO	Cancer association	Cancer	Ref
HID1 (DMC1)	139	chr17q25.1	Transmembrane	Downregulated	Breast, cervix, liver, lung,	[53]
					thyroid, stomach, kidney
FER1L4	135	20q11.22	lncRNA	Downregulated	Stomach	[54,81]
TTLL3	128	3p25.3	Tubulin-tyrosine	Downregulated	Colon	[55,82]
			ligase activity
RIF1	126	2q23.3	DNA repair	Anti-apoptotic	Breast	[56]
KLHDC7A (FLJ38753)	116	1p36.13	Transmembrane	-	-	-
SBNO1	110	12q24.31	DNA binding	Proliferation	Lung	[57]
B) Bead UC GRN
Gene	Degree	Locus	Description/GO	Cancer association	Cancer	Ref
MGC15885	119	15q22.2	ncRNA	-	-	-
RNF17 (TDRD4)	107	13q12.12	Spermatid development	Potential cancer CT antigen	Liver	[58]
OR6S1	105	14q11.2	G-protein coupled	-	-	-
			receptor signaling pathway
ACTR3BP5	104	10p11.1	Pseudogene	-	-	-
TPTE2P3	98	13q14.3	Pseudogene	-	-	-
TMED7	90	5q22.3	Protein transport	Upregulated	Nasopharynx carcinoma	[59]
					cell line
C) Oligo UC GRN
Gene	Degree	Locus	Description/GO	Cancer association	Cancer	Ref
CYP4A11	184	1p33	Monooxygenase activity	Promotes angiogenesis	Lung	[60]
				and metastasis
GJC2	162	1q42.13	Gap junction channel	-	-	-
			activity
GPATCH4	147	1q23.1	Nucleic acid binding	-	-	-
ADAM5	141	8p11.22	Metalloendopeptidase	-	-	-
			activity, pseudogene
DKFZP434A062	120	9q34.3	Uncharacterized protein	-	-	-
SLC38A3 (SNAT3)	118	3p21.31	Symporter activity	malignancy marker	Glioma	[61]

Shown are the gene symbols of the hub genes of the (A) RNAseq, (B) Bead and (C) Oligo UC GRN, their number of interactors (Degree), chromosomal location (Locus), functional description when available from GO or gene description, literature-based evidence or property for a cancer association, cancer types (Cancer) and literature citation (Ref) when available.

Overlapped only for a single term, i.e., for the regulation of neuron differentiation (GO:0045664). For the RNAseq and Oligo GRN we observed nine terms in agreement, e.g., G-protein coupled receptor signaling pathway (GO:0007186), ion transport (GO:0006811) and sensory perception (GO:0007600). We observed that the average degree centrality of the networks with randomized gene labels was similar across the gene sets and independent of the number of genes of a gene set. However, our RNAseq dataset mostly considered muscle-invasive UCs and shows terms associated to invasiveness. The Bead data comprised terms that were nuclear while the terms predicted from the Oligo data were predominantly extracellular membrane associated.

In this section we quantified the extend of local and globally shared interactions between the inferred GRNs and ppi networks. We assessed the global extent of interactions from PPI databases that were present in the UC GRNs by comparing the entire GRNs networks to String, Biogrid, Hprd, Intact, Mint, Graphite, CPDB and SingNet (Additional file 1: Table S10).

As observed on the subnetwork level using all String interactions the F-scores between the RNAseq and the Oligo UC GRN were similar. In contrast the Bead UC GRN showed a lower F-score for all PPI databases compared to the RNA-seq GRN and Bead GRN (Additional file 1: Table S10).

Further, we compared the relative quantity of PPI interactions for the identified significant subnetworks of the gene sets for biological processes, genomic co-located genes and gene families. To avoid the comparison to subnetworks with no known protein-protein associations we used String as reference as it was the largest collection of PPI interactions that we considered in our study. For each gene set we computed F-scores by comparing the corresponding subnetwork of the GRN to the subnetwork of the String network reference. Figure 4 shows the cumulative F-score distributions between the RNAseq, Bead and Oligo GRN separately for the 299 GO terms, 28 gene family subnetworks and 40 chromosomal 1 Mb regions. In addition, we repeated the analysis for each GRN 25 times using a reference subnetwork where the gene labels were randomized (Figure 4).

Figure 4

Network inference performance using interactions from the STRING database. Shown are the F-score distributions for commonly significant GRN subnetworks of the GPEA analysis for A) 299 Gene Ontology Biological Process terms; B) co-located genes of 1 Mb genomic regions; C) 28 gene families; and D) p-value distribution (FDR) estimated for all gene sets.

The average F-scores were highest for the subnetworks of genomic co-located genes (F=0.26, F=0.18 and F=0.26) and for the gene family gene sets (F=0.15, F=0.12, F=0.15). The Gene Ontology gene sets had the lowest average F-scores compared to the genomic and gene family subnetworks (F=0.044, F=0.029, F=0.039) (Figure 4 and Additional file 1: Table S11). The observations are in agreement with the global analysis for RNAseq and Oligo GRN, where the Bead UC GRN has the tendency to perform worse. However, the RNAseq UC GRN shows the tendency for smaller p-values over the Bead and Oligo UC GRN (Figure 4D).

The reference network with randomized gene labels were significantly lower compared to the GRN for all comparisons. The Bead GRN had the tendency to show a significantly lower mean F-score compared to the RNAseq and Oligo GRN for the Gene Ontology (t-test, p=0.000007, p=0.002390) and the genomic 1Mb window subnetworks (p=0.000675, p=0.001697). For the gene family subnetworks we did not observe a significant difference among the three GRN.

Discussion

In this paper, we have presented novel perspectives and applications for the identification of UC molecular targets using GRNs. Specifically, we performed a structural, functional and comparative analysis across three UC GRNs that were independently inferred from three large-scale RNAseq, Bead and Oligo gene expression datasets. Our objective was to identify putative prognostic UC targets for a subsequent investigation in UC tumors on the basis of their enrichment in functional subnetworks and hub gene analysis. Our results demonstrate that GRNs are highly dataset specific on the gene interaction level and showed large similarities across the functional subnetwork levels. The RNAseq based GRN showed the most prominent functional enrichment and is thus the data type of choice for a network inference. The RNAseq and Oligo GRN showed a similar inference performance based on public interaction databases and outperformed the Bead based GRN.

On the structural level, the three inferred GRNs were observed to follow a power law distribution [62] that is common for inferred and experimental biological networks [63-65]. Our results demonstrated that the network structure at the gene level of GRNs are highly dependent on the individual gene expression dataset. On the gene interaction level the pairwise comparison between the networks showed only an overlap of 2% and only 0.25% of all interactions are common among the three networks (Figure 1). There are three main explanations for this observation. The first reason is that the BC3Net algorithm considers only the strongest interaction neighbors for each gene and is thus highly dependent on the search space of the genes that are included in the dataset. The second reason is that the variations caused by concordance differences of the expression are dependent on technical properties of the individual gene expression platforms and platform dependent data processing procedures. The third reason is that the datasets represent varying proportions of different tumor grades and stages from individual patients that represent a complex condition phenotype. Further, gene expression profiles of tumor tissues are highly heterogeneous on the molecular and tissue level, i.e., tumor clonal variation within and between different patient samples [66].

In [67] a guilt-by-association approach was developped to predict molecular roles of genes with unknown functions. The “guilt-by-association” property of genes that are connected in a defined network can also be used for a functional enrichment analysis for gene pairs which have known functions and are involved in the same biological processes. We identified significant functional GRN subnetworks by performing a gene pair enrichment analysis (GPEA) for defined gene sets. We used the terminology gene pair enrichment analysis to distinguish the latter from the terminology for a gene-based enrichment analysis which has no structural constraint. The concept for the analysis was introduced from graph theory [68] and has been developped and applied for the identification of significant protein complex and ontology gene sets in PPI and inferred networks [69-72].

A total of 5 to 10% of all tested Gene Ontology Biological process terms, 2 to 10% of gene sets of co-located genes and 9 to 20% gene family gene sets showed a significant subnetworks by the enrichment of inferred interactions (Figure 1). RNAseq based network showed more than twice the proportion of significant subnetworks compared to the Bead and Oligo microarray based GRN. Our results showed in a quantative manner that RNAseq is beneficial for GRN inference compared to Bead and Oligo microarray based data. The major advantages of RNAseq are more accurate measurement of the dynamic range of low and highly expressed genes [73] and thus gives a better resolution of the underlying functional dependency structure between the genes.

In contrast to the low similarity that was observed between the GRNs on the structural interaction level, we observed high similarities on the functional subnetwork level (Figure 1). The fraction of significant biological process Gene Ontology terms that were common across the three UC GRN was above 55%. For the gene family subnetworks we observed a similarity for 25.45% and the lowest percentage for genomic co-located gene subnetworks 6.8% (Figure 1).

The networks described a prominent enrichment for known cancer hallmarks [74] with significant GO subnetworks related to immune response, cell cycle, signal transduction, DNA repair, translation, proteolysis, metabolic terms such as respiration and cell morphogenesis, adhesion and migration. Further, over 50% of the significant GO subnetworks were highly enriched by known cancer genes defined by the cosmic cancer census [31] across the three UC GRNs. We observed that the fraction of cancer associated subnetworks is prominently lower in relevance network inference methods. This may result from low dependency measures of relevant interactions of genes in a more complex context being excluded from a relevance network by a global threshold. For other GRN inference methods we expect similar results to the results presented by the BC3Net that is based on the C3Net. A C3Net infers a core structure of a GRN and thus infers only a subnetwork of other GRN inference methods based on mutual information [10]. For each gene in a C3Net at most one gene neighbor with strongest mutual dependency is selected, which results in a highly reduced time complexity for multiple testing of mutual information. The C3Net and BC3Net GRNs inference method is therefore less time consuming which makes the inference of very large GRNs (>20K genes) feasible in a reasonable time.

On the genomic level, the GRNs were investigated for genomic UC targets, where we identified genomic regions with known diagnostic and prognostic properties for urothelial cancer such as 1q23.3 [52] (Oligo GRN), 8q24.3 [51] (Bead GRN) and 5q31.3 [48-50] (RNAseq GRN). The identified genomic regions can link to chromosomal aberrations, histone modifications, changes in epigenetic regulation (methylation), regulatory elements and spatial chromosome organization in the nucleus. These processes are commonly deregulated in cancer. For example the impairment of DNA repair mechanisms leads to an accumulation of chromosomal aberrations that are frequently observed in the progression of UC [75]. The identification of subnetworks of genomic regions from co-located genes therefore provided a powerful tool to identify putative novel genomic targets from cancer gene expression datasets.

In the analysis of gene families we found CD molecules as the most prominent gene family. CD molecules are promising targets for novel cancer immunotherapies such as CD47 [76]. Some popular UC biomarkers target proteins of an entire gene family and not a single gene product such as Keratins [77] and Kalikrein proteins [78] which are popular tumor markers for UC. Gene families are crucial in cancer research [79] because they represent groups of genes that are functionally highly redundant and represent potential targets of the underlying molecular heterogeneity that is observed for malignant processes. We showed that the identification and ranking of functional and co-located gene sets and gene families using our GPEA on GRNs is a versatile approach for the generation of novel targets and molecular understanding of the properties of urothelial cancer from the perspective of large-scale tumor tissue gene expression data.

Hub genes of GRNs reflect the most prominent dependencies of the expression profile to a large number of genes. We identified hub genes such as HID1 [53] (RNAseq) RNF17 (TDRD4) [58] (Bead), CYP4A11 [60] (Oligo) for the individual GRNs which show in the literature strong evidence for cancer related diagnostic and prognostic properties. Further, we performed a degree centrality analysis of the GRNs that showed that the degree centrality of the genes allow to target promising mediators of cancer related cellular activities and signaling processes.

In addition, we performed a quantitative comparison of protein interactions for the RNAseq, Bead and Oligo UC GRNs. We note that the overlap for protein interaction data and GRNs is expected to be low and non-random. For example the most prominent PPI interactions that can be found in a GRNs are physical interactions of genes corresponding to large protein complexes (e.g. ribosome biogenesis and proteasome) in contrast to more transient protein interactions [11,80]. A GRN is inferred from gene expression data and thus can only detect indirect association to the protein level of a gene network. However, the analysis allowed to compare network properties between the Oligo, Bead and RNAseq data and pointed to the tendency that Oligo expression data should be prefered over Bead expression data for a GRN inference.

Conclusion

On the functional and structural level our results demonstrated that RNAseq based data is the preferred data type for a GRN inference. GRNs are highly dataset-specific at the interaction level, while at the global functional level they are highly similar. GRN inference is a powerful tool to provide a database of novel UC targets that can be studied for prognostic and diagnostic clinical applications [48,49,52,58,60].

Supplementary Materials: Urothelial cancer gene regulatory networks inferred from large-scale RNAseq, Bead and Oligo gene expression data.