Analysis of molecular intra-patient variation and delineation of a prognostic 12-gene signature in non-muscle invasive bladder cancer; technology transfer from microarrays to PCR.

BACKGROUND: Multiple clinical risk factors and genetic profiles have been demonstrated to predict progression of non-muscle invasive bladder cancer; however, no easily clinical applicable gene signature has been developed to predict disease progression independent of disease stage and grade.
METHODS: We measured the intra-patient variation of an 88-gene progression signature using 39 metachronous tumours from 17 patients. For delineation of the optimal quantitative reverse transcriptase PCR panel of markers, we used 115 tumour samples from patients in Denmark, Sweden, UK and Spain.
RESULTS: Analysis of intra-patient variation of the molecular markers showed 71% similar classification results. A final panel of 12 genes was selected, showing significant correlation with outcome. In multivariate Cox regression analysis, we found that the 12-gene signature was an independent prognostic factor (hazard ratio=7.4 (95% confidence interval: 3.4-15.9), P<0.001) when adjusting for stage, grade and treatment. Independent validation of the 12-gene panel and the determined cut-off values is needed and ongoing.
CONCLUSION: Intra-patient marker variation in metachronous tumours is present. Therefore, to increase test sensitivity, it may be necessary to test several metachronous tumours from a patient's disease course. A PCR-based 12-gene signature significantly predicts disease progression in patients with non-muscle invasive bladder cancer.

A total of 70 530 new bladder cancer cases and 14 680 deaths were estimated in the United States alone in 2010 (Jemal ). About 75% of patients present with non-muscle invasive tumours were mostly treated with a local, organ-sparing approach. The remaining 25% of patients are initially diagnosed with muscle invasive cancers, most often requiring immediate and more radical treatment. More than 60% of patients with non-muscle invasive bladder cancer experience tumour recurrences and around 20% of the patients develop disease progression to a muscle invasive bladder cancer (Millan-Rodriguez ; Sylvester ). Clinical risk factors associated with a high risk of disease progression to a muscle invasive cancer include deep invasion of the lamina propria, high-grade tumour, large tumour size, concurrent carcinoma in situ (CIS), tumour multiplicity and recurrence of high-risk non-muscle invasive tumours (Hermann ). The clinical risk factors cannot predict the individual disease course accurately, and currently no molecular markers are available for clinical use for predicting tumour recurrence or later disease progression (Ehdaie and Theodorescu, 2008).

Several microarray-based gene-expression signatures for diagnosis or prognosis have been identified in many cancers, including bladder cancer (Rosenwald ; van de Vijver ; van’t Veer ; Ramaswamy ; Dyrskjot ). Although clinically applicable, it may, however, be easier to implement tests for clinical use when simpler, cheaper and more sensitive methods are being used. We have previously identified a molecular gene-expression signature for predicting disease progression in patients with non-muscle invasive bladder cancer (Dyrskjot ). The progression signature was validated successfully in a retrospective study that included 294 patients from five different countries. The molecular signature was significantly correlated with progression-free survival (P<0.001) and it was shown to be significantly associated with disease progression after adjustment for age, sex, stage, grade and treatment (hazard ratio (HR)=2.3; P=0.007) (Dyrskjot ). One of the objectives to be addressed before clinical implementation is intra-patient variation of the molecular signature when applied to multiple tumours from the same patient. How often do we need to test in each patient? Another objective is the generation of a clinically applicable gene signature, and this requires transferring the markers from the microarray platform to a robust PCR-based platform.

Here we address the intra-patient variation of the original 88-gene signature by analysing 39 metachronous tumours from 17 patients. We combined this analysis with mutation screening of the genes frequently mutated in bladder cancer: FGFR3 (16–74% mutated – dependent on stage), PIK3CA (13–27% mutated) and RAS family (NRAS, KRAS and HRAS; 13% mutated) (Billerey ; Kompier ). Then, we describe the transfer of the optimal signature genes from the microarray platform, plus additional prognostic genes identified, to a PCR platform for easier clinical implementation of the molecular classifier.

Materials and methods

Patients and tumour samples

The majority of the patients and tumour samples used in this work were used previously for validating our microarray-based 88-gene progression signature (Dyrskjot ). Urothelial carcinomas were taken from patients that were operated in the years 1987–2000 in hospitals in Denmark, Sweden, Spain and England and frozen immediately at −80 °C. All patients were followed according to the valid guidelines at patient inclusion time. Informed written consent was obtained from all patients, and research protocols were approved by the institutional review boards or ethical committees in all involved countries. Diagnostic pathology slides were reevaluated according to the World Health Organisation (WHO) 2004 guidelines. However, for 19 out of the 115 patients, it was not possible to acquire the diagnostic sections for review and consequently the original grading was translated into the WHO 2004 grading system (G1+G2=low grade, G3+G4=high grade), although this translation is too simplified in some cases. Progression-free survival time was recorded from sampling visit and censored at the time of the last control cystoscopy or at cystectomy. Progression of the disease was defined as (i) histologically confirmed invasion into the bladder muscle or (ii) cases where the patient died of bladder cancer without a verified progression event as described above (n=5). In case of bladder cancer death, survival was recorded from the date at which tissue was collected and until the last annotation of the patient being alive.

RNA extraction and cDNA generation

Total RNA was extracted from the Danish, and English samples using a standard Trizol RNA extraction method (Invitrogen, Carlsbad, CA, USA). Total RNA from Swedish and Spanish samples was extracted using RNeasy mini kit (Qiagen, Valencia, CA, USA). All RNA was quality controlled using an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) (criteria: 28S/18S >1 and RIN >5). In all, 1 μg of total RNA was DNAse-treated using DNase I Amplification Grade (Invitrogen) to avoid amplification of genomic DNA. Adequate DNAse treatment was verified in all samples by quantitative reverse transcriptase PCR (qRT–PCR) using intron spanning GAPDH primers. The DNAse-treated total RNA was converted to cDNA using oligo (dT) priming and SuperScript II Reverse Transcriptase (Invitrogen) according to the standard protocols.

Initial selection of 96 genes for analysis

Genes were selected using two different approaches. First, genes were selected from the prognostic gene signatures published previously (Dyrskjot ). Second, we reanalysed previously published microarray data (Dyrskjot , 2004; Aaboe ; Dyrskjot ). For each gene, the microarray hybridisation intensity was plotted against progression-free survival of each patient and linear and quadratic regression equations were computed. The noise was calculated as a function of squared deviations of each patient from the squared regression line. Then, the ratio of linear regression coefficient to noise was calculated for each gene and used as a measure of ‘signal-to-noise’ ratio. Genes were then selected based on ‘signal-to-noise’ performance. The two different approaches showed large overlap in genes selected, and we selected those showing the best performance with each approach.

Primer design and PCR assays

Microarray layout files were used to map the Affymetrix (Santa Clara, CA, USA) probe sets to mRNA sequences and this region was then used to design primers for qRT–PCR assays. For each of the selected candidate genes, we designed three primer pairs using PRIMER3 software (Rozen and Skaletsky, 2000). For sensitivity measures, we used 1, 0.1 and 0.01 ng cDNA template. A minimum of two replicates had to amplify at a particular concentration to qualify as the limit of sensitivity. Furthermore, the threshold cycle (Ct) spread between the replicates at this concentration had to be less than one cycle. Specificity measures were performed using no template control measurements where Ct should be >40. The qRT–PCR using the final primer set was carried out using a 7900HT Fast Real-Time PCR System (Life Technologies, Carlsbad, CA, USA) in 384-well plates. A total of 115 tumour samples were analysed. All reactions were performed in triplicate in 10-μl volumes using SYBR Green PCR Master Mix (Life Technologies) and the standard protocols were applied.

Normalisation, optimal 12-gene set selection and molecular classification

To identify the best-performing candidate genes, we normalised data using average total-patient Ct values and using an average global total. Classifier genes were selected based on different statistical methods applied to normalised Ct values and clinical outcome (progression vs no progression): Pearson’s correlation, P<0.01; Wilcoxon signed-rank test, P<0.01; Cox regression analysis, P<0.01; and area under curve (AUC) from receiver operating characteristics (ROCs), AUC>0.65. Classifier genes with high correlation were not excluded from the final signature. Computation of risk scores utilised non-normalised Ct values based on the formula: average (Ct (genes downregulated in progressing tumours)−average (Ct (genes upregulated in progressing tumours)). This approach eliminated the requirement of using normalised Ct values. The optimal cut-off value (risk score with the highest accuracy) and the 90% sensitivity cut-off value were identified from ROC curves, and risk scores were then dichotomised based on these two different cut-off values. Dichotomised signature values were used to classify samples as either high or low risk of progression (1=high risk and 0=low risk).

Affymetrix Exon ST 1.0 array analysis and classifier construction

Affymetrix human Exon ST 1.0 arrays were used for measuring gene expression in 39 tumours from 17 patients as described previously (Thorsen ). itPLIER normalisation and generation of gene-expression measures were performed in GeneSpring GX 10.0 with no baseline adjustment. To identify the original 88-gene signature, we used BLASTN (NCBI) against human RefSeq RNA sequences for all original probe sequences to updated Gene Symbol annotation. Following, we used the Gene Symbols for identifying the corresponding probe sets on the Exon ST 1.0 arrays. We preferentially selected ‘core’ probe set with the highest number of exon clusters. A maximum likelihood classifier was constructed as described previously (Dyrskjot ) using 10 tumours that were originally used for building the classifier and for testing the 39 independent tumour samples.

Mutation analysis

Mutation analysis of the genes FGFR3, PIK3CA, NRAS, KRAS and HRAS was performed essentially as described previously (Kompier ) for 36 out of the 39 tumours samples analysed on the Exon ST 1.0 microarray platform.

Statistical procedures

We used STATA 10.0 (Stata Corporation, College Station, TX, USA) for statistical analysis. Variables with a P-value <0.05 in univariate analysis were included in multivariate analysis. The assumptions of proportional hazards were verified.

Results

Intra-patient variation

We set out to analyse the progression signature in 39 metachronous bladder tumours from 17 patients by Exon ST 1.0 microarray analysis to measure the intra-patient variation of the signature. The mutational status of the samples was also analysed to obtain an indication of the genetic relatedness of the metachronous tumours. The classification result and mutation analysis of the genes FGFR3, PIK3CA, NRAS, KRAS and HRAS are shown in Table 1. NRAS, HRAS and KRAS were not mutated in any of the tumours analysed. Multiple tumours from 12 patients (71%) showed similar intra-patient classification results, however, in another five cases the intra-patient analysis results were dissimilar; two of these patients showed variation in mutation status of the FGFR3 gene, indicating differences in tumour clonality. All other tumours showed similar mutation status.

Table 1

Intra-patient variation based on 88-gene classifier results and mutation analysis

				Mutation analysis
Sample	Stage	Grade a	88-gene classifier b	FGFR3	PIK3CA	NRAS	HRAS	KRAS
532-12	Ta	II	1	S249C	WT	WT	WT	WT
532-15	Ta	II	1	S249C	WT	WT	WT	WT
1888-1	Ta	II	1	S249C	WT	WT	WT	WT
1888-2	Ta	II	1	S249C	WT	WT	WT	WT
1841-4	Ta	II	1	S249C	WT	WT	WT	WT
1841-6	Ta	II	1	S249C	WT	WT	WT	WT
1926-1	Ta	II	1	G372C	WT	WT	WT	WT
1926-4	Ta	II	1	G372C	WT	WT	WT	WT
1926-6	Ta	II	1	G372C	WT	WT	WT	WT
1657-1	Ta	I	1	S249C	WT	WT	WT	WT
1657-3	Ta	II	1	S249C	WT	WT	WT	WT
1899-2	Ta	III	1	G372C	WT	WT	WT	WT
1899-3	T2–4	III	2c	WT	WT	WT	WT	WT
1992-1	Ta	I	2	S249C	WT	WT	WT	WT
1992-3	T1	II	1	S249C	WT	WT	WT	WT
2010-1	T1	III	2	S249C	WT	WT	WT	WT
2010-4	T2–4	III	2	S249C	WT	WT	WT	WT
1375-1	T1	II	1	S249C	E545K	WT	WT	WT
1375-7	Ta	II	1	S249C	E545K	WT	WT	WT
1735-1	T1	III	2	WT	WT	WT	WT	WT
1735-4	Ta	II	2	WT	WT	WT	WT	WT
1702-1	Ta	I	1	S249C	WT	WT	WT	WT
1702-7	Ta	II	1	S249C	WT	WT	WT	WT
1702-11	Ta	I	1	S249C	WT	WT	WT	WT
1790-1	Ta	I	1	Y375C	WT	WT	WT	WT
1790-3	Ta	II	1	Y375C	WT	WT	WT	WT
562-3	Ta	I	2	Y375C	WT	WT	WT	WT
562-5	Ta	II	2	Y375C	WT	WT	WT	WT
1521-1	Ta	I	1	S249C	WT	WT	WT	WT
1521-6	Ta	I	1	S249C	WT	WT	WT	WT
1521-10	T1	II	1	—	—	—	—	—
1732-1	T1	II	2	—	—	—	—	—
1732-11	Ta	II	1	S249C	WT	WT	WT	WT
1732-12	T2–4	III	2	—	—	—	—	—
1829-1	Ta	II	1	S249C	WT	WT	WT	WT
1829-2	Ta	II	2	S249C	WT	WT	WT	WT
1829-3	T1	III	1	S249C	WT	WT	WT	WT
2363-1	T1	III	1	S249C	WT	WT	WT	WT
2363-2	T2–4	III	2	WT	WT	WT	WT	WT

Bergkvist grading.

1=classified as a non-progressing tumour and 2=classified as a progressing tumour.

Less than 5% difference in classification distance between the two classes. Bold denotes intra-patient differences in 88-gene classifier results or in mutation status.

Patient and marker selection

Initially, we selected 96 genes associated with disease aggressiveness from published gene signatures (38 genes) and from reanalysis of microarray data (58 genes). Following primer design and initial validation of primer sensitivity and specificity, we decreased the number of interesting genes to 35 in total (Table 2). We performed qRT–PCR measurements of the 35 genes using 115 tumour samples from patients in Denmark, Sweden, UK and Spain. Previously, 102 of these tumours were used for validating our microarray-based 88-gene signature. Patients with no disease progression were followed for a minimum of 60 months (the median follow-up time was 88 months). Clinical and histopathological characteristics for all patients are shown in Table 3.

Table 2

Primer sequences for the 35 candidate genes

Gene target a	Upstream primer sequence 5′–3′	Downstream primer sequence 5′–3′	Gene selection b
ADAM10	5′-CAGTATTACTTATGGGAATTGCTCTG-3′	5′-TTGGATTACTACTTGGAGTATGAACAC-3′	1
AURKB	5′-TTTGCTATGAGCTGCTGGTG-3′	5′-ACTTTAGGTCCACCTTGACGATG-3′	1
BIRC5	5′-CTGAAGTCTGGCGTAAGATGATG-3′	5′-GAAGCTGTAACAATCCACCCTG-3′	1
C10orf58	5′-GTAAACCTACTTTCTGTTCTGGAAGC-3′	5′-TTTTCTCTGAGGCCAAAGTCTG-3′	1
CDC20	5′-AGTCCAATGTCCTGGCAACAG-3′	5′-CCAGAGCACACATTCCAGATG-3′	1
CDC25B	5′-GATGGAAGGTTGGATGGATG-3′	5′-ACCTGGTTTGGGTATGCAAG-3′	1
COL4A1	5′-CTGCCTGGAGGAGTTTAGAAGTG-3′	5′-CTGTAAGCGTTTGCGTAGTAATTG-3′	1
FABP4	5′-AGAGAAAACGAGAGGATGATAAACTG-3′	5′-CTTATGCTCTCTCATAAACTCTCGTG-3′	1
IGF2	5′-CATCGTTGAGGAGTGCTGTTTC-3′	5′-GGGTAGCACAGTACGTCTCCAG-3′	1
ITGB4	5′-CATCATCCCTGACATCCCTATC-3′	5′-GTAGAACGTCATCGCTGTACATAAG-3′	1
KPNA2	5′-GCAGATTTTAAGACACAAAAGGAAG-3′	5′-AAGGTACACAATCTGTTCAACTGTTC-3′	1
MAT2B	5′-TTGTCTAAAGAAACTAAAGGGCAGTC-3′	5′-AGTTTAGCCAGGACAAACAAAATG-3′	1
MBNL2	5′-ACTTCATCCAGTGCCCACTTTC-3′	5′-GGGGTTACAGGTGCTAGGTAAGG-3′	1
MCM7	5′-GAGATGTCAAAGGACTCTCTTCTAGG-3′	5′-GCAAATATCACATCTGCTGGTC-3′	1
MSN	5′-CCTGACCTTGAGGAGTCTTGTG-3′	5′-AATATAGGACATATCACCAAGTGAGC-3′	1
PPP2R5C	5′-GTACTACATTGAAAATAAACCGGTGAC-3′	5′-TACATTTTGGAAAGAGTGAAGATGC-3′	1
TCF4	5′-GAATCACATGGGACAGATGTAAAAG-3′	5′-AATACAGCTGTTAAGGAAGTGGTCTC-3′	1
ACTA2	5′-GTCTCTAGCACACAACTGTGAATGTC-3′	5′-CTAGGAATGATTTGGAAAAGAACTG-3′	2
CDH5	5′-AAACAATTCCTGTAACCTTCTATTTTC-3′	5′-CTTGTCATGCACCAGTTTGG-3′	2
CDKN3	5′-ATCTCTACCAGCAATGTGGAATTATC-3′	5′-CTATGTCAGGAGTCCCTCCATC-3′	2
COL18A1	5′-GGGCTGGTTCTGTAATTGTGTG-3′	5′-AAAAGGTCACTTTTATTTGCCTGTC-3′	2
COL4A3BP	5′-TTTCTGTGGATCATGACAGTGC-3′	5′-CAAGGTTTGACAAATCATAGCAAC-3′	2
CPS1	5′-GGAAGTAAGGTTCATTCCCTTAAGAC-3′	5′-CCTTACAGTGGGTGGAATTAATAGTG-3′	2
DCTD	5′-TGTGCCCCTTCTCTTTAATCTC-3′	5′-GAAAGCCTTTTCTCAACACAGG-3′	2
IER2	5′-CTTGCCAGGGAGTTTCTGAG-3′	5′-ATTTCTAACAAAACGCCAGGTAGAC-3′	2
LBR	5′-ATCAGAAAGTGGTGGCGTTTTC-3′	5′-TTACCAGGGAAAGAATTTAATGTCC-3′	2
LGALS1	5′-CTGAATCTCAAACCTGGAGAGTG-3′	5′-GGTTCAGCACGAAGCTCTTAGC-3′	2
NEK1	5′-CTAAAAGACCAGCTTCAGGACAAAAC-3′	5′-CTAAAGGTATTCCATATTTAGCGGC-3′	2
NR1H3	5′-GGAATTCATCAACCCCATCTTC-3′	5′-GATAGCAATGAGCAAGGCAAAC-3′	2
PEA15	5′-ACTCCTTATATTGCTGTGAGATTGC-3′	5′-ACCTTTATTCCGGGTTAGAACAAG-3′	2
SKAP2	5′-TGGAGATGTATGATATTTGAGAGTCC-3′	5′-CTAAATCCAAAGCATTTGCAGAC-3′	2
SDC1	5′-AGACACCTTGGACATCCTCCTC-3′	5′-TAAGCAAGTAAGTGCAGGAGCC-3′	2
SEC14L1	5′-TGTTTCTACCTTTAGTACCTTGCCAC-3′	5′-AGTACTAAGAAATGGGAAATGACAGC-3′	2
UBE2C	5′-TCTAGGAGAACCCAACATTGATAGTC-3′	5′-TCTTGCAGGTACTTCTTAAAAGCTG-3′	2
WNT2B	5′-ATAAGAAACTGTGCAAGCTCCC-3′	5′-TCTACTCTCCCTTCAAATCTCCAG-3′	2

Genes in bold are included in the final 12-gene signature.

Gene selection method: 1=from previously published gene signatures and 2=from reanalysis of gene-expression data.

Table 3

Clinical and histopathological characteristics for all 115 patients

Number of patients	115
Median follow-up time in months for all patients (range)	69 (1–216)
Median follow-up time in months for progressing patients (range)	20 (1–123)
Median follow-up time in months for non-progressing patients (range)	88 (61–216)
Median age (range)	70 (33–88)
Male–female ratio	4.9
Stage
Ta	62 (54%)
T1	53 (46%)
Grading (WHO 2004)
PUNLMP	17 (15%)
Low grade	33 (29%)
High grade	65 (56%)
Concomitant CIS
Yes	19 (17%)
No	50 (43%)
Unspecified	46 (40%)
Adjuvant therapy (BCG or MMC)
Yes	34 (30%)
No	81 (70%)
Number of progression events to stage T2–4 bladder cancer
Ta	13 (30%)
T1	33 (70%)

Abbreviations: BCG=bacillus Calmette-Guerin; CIS=concurrent carcinoma in situ; MMC=mitomycin C; WHO=World Health Organisation.

12-gene PCR signature

We selected the genes for inclusion in the optimal PCR signature by comparing normalised Ct values and clinical outcome using Pearson’s correlations, ROC analysis, Wilcoxon signed-rank tests and Cox regression analysis. Using these different statistical methods in a combined approach, we selected 12 genes that showed significant correlation with outcome (see Figure 1A and Table 4); seven were upregulated in progressing tumours compared with non-progressing tumours (KPNA2, BIRC5, UBE2C*, CDC25B, COL4A1*, MSN and COL18A1*) and five were downregulated in progressing tumours compared with non-progressing tumours (COL4A3BP*, MBNL2, NEK1*, FABP4 and SKAP2*). The six genes marked with an asterisk were not previously included in our gene signature for outcome prediction. Selecting fewer genes using more stringent selection criteria did not result in increased classification performance (results not shown). A schematic representation of the steps involved in reaching the final 12-gene PCR signature is shown in Figure 2.

Figure 1

The 12-gene PCR signature for predicting disease progression. (A) Expression patterns of the 12 genes in the 115 tumour samples analysed. The upper seven genes are upregulated in progressing tumours and the lower five genes are downregulated in progressing tumours. Yellow, upregulation of the gene; blue, downregulation; and black, median expression. The vertical red line separates tumours classified with the progression signature (right) from tumours with the non-progression signature (left) based on the optimal cut-off value. The black boxes below the heat map denote progressing tumours, stage T1 tumours and high-grade tumours. (B) ROC curve for the 12-gene signature values (average Ct (genes downregulated in progressing tumours)−average Ct (genes upregulated in progressing tumours)). The value with the highest accuracy is indicated as the optimal cut-off value for the dichotomised test (used in C, D and E), and the 90% sensitivity value is indicated (used in F, G and H). (C, F) Kaplan–Meier plots of progression-free survival as a function of the 12-gene PCR signature. All 115 patients are included in the analysis. (D, G) Kaplan–Meier plots of progression-free survival as a function of the 12-gene PCR signature when including high-risk patients only (stage T1 or high grade or CIS). (E, H) Kaplan–Meier plots of progression-free survival as a function of the 12-gene PCR signature when including low-risk patients only.

Table 4

Progression-associated genes selected for the final 12-gene PCR signature and associated statistics based on 115 tumour samples

			Cox regression analysis
GENE	Pearson’s correlation a	Wilcoxon signed-rank test P -value	Beta coefficient	P -value	ROC AUC
CDC25B	−0.268	0.004	−0.605	0.002	0.665
KPNA2	−0.234	0.004	−0.490	0.006	0.663
BIRC5	−0.251	0.004	−0.449	0.005	0.667
COL18A1	−0.307	<0.001	−0.368	0.002	0.704
MSN	−0.277	<0.001	−0.328	0.007	0.690
UBE2C	−0.331	<0.001	−0.315	<0.001	0.729
COL4A1	−0.345	<0.001	−0.300	<0.001	0.740
FABP4	0.374	<0.001	0.193	<0.001	0.729
MBNL2	0.297	<0.001	0.459	0.001	0.707
SKAP2	0.270	0.004	0.566	0.004	0.664
COL4A3BP	0.359	<0.001	0.777	<0.001	0.746
NEK1	0.253	0.003	0.804	0.003	0.671

Abbreviations: AUC=area under curve; ROC=receiver operating characteristic.

P<0.01 if the Pearson’s correlation is >0.176 or <−0.176. A negative value denotes genes upregulated in progressing tumours.

Figure 2

Schematic overview of steps involved in delineation of the final 12-gene PCR signature.

Molecular classification

Molecular prediction of outcome was carried out for all 115 patients based on the generated PCR risk scores. An optimal cut-off value (0.79; 76% sensitivity and 86% specificity) was identified using ROC analysis (Figure 1B) and the dichotomised scores showed significant correlation with progression-free survival (P<0.0001, log-rank test; Figure 1C). In multivariate Cox regression analysis, we found that the 12-gene signature was an independent prognostic variable (HR=7.4 (95% confidence interval: 3.4–15.9), P<0.001) when adjusting for stage, grade and BCG/MMC treatment (Table 5). When only including high-risk patients in the analysis (i.e., stage T1, high grade or concomitant CIS), we also found a significant correlation with progression-free survival (P<0.0001, log-rank test; Figure 1D). Similarly, when including only the low-risk subpopulation in the analysis, we also found a significant correlation with progression-free survival (P=0.0002, log-rank test; Figure 1E). We calculated the predictive accuracy using Harrell’s C to 73% for clinical variables only, 75% for the 12-gene signature, and this number was raised to 82% when including clinical variables and the 12-gene signature in the model. Consequently, accurate prediction of progression is increased when including the molecular markers. When changing the cut-off value to −0.17 to ensure a 90% sensitivity (and 61% specificity) of the test, we also observed significant correlations with outcome for the entire patient cohort and the risk subgroups as shown in Figure 1F–H.

Table 5

Univariate and multivariate Cox regression analysis of progression-free survival

	Univariate analysis		Multivariate analysis
	Hazard ratio (95% CI)	P -value	Hazard ratio (95% CI)	P -value
Age (per 5 year increment)	1.21 (1.02–1.43)	0.024	1.02 (0.85–1.23)	0.806
Sex (men vs women)	0.62 (0.31–1.25)	0.184
Stage (T1 vs Ta)	4.07 (2.13–7.77)	<0.001	1.64 (0.63–4.28)	0.311
Grade (high vs low+PUNLMP)	4.82 (2.25–10.34)	<0.001	1.25 (0.41–3.86)	0.703
Concomitant CISa (presence vs absence)	2.37 (1.06–5.30)	0.035
Treatment (BCG and MMC vs no treatment)	0.23 (0.09–0.58)	0.002	0.20 (0.07–0.57)	0.003
12-gene PCR signature (high risk vs low risk)	9.36 (4.66–18.81)	<0.001	8.30 (3.66–18.84)	<0.001

Abbreviations: BCG=bacillus Calmette-Guerin; CI=confidence interval; CIS=concurrent carcinoma in situ; MMC=mitomycin C.

CIS status was only known for 69 patients and consequently not included in the multivariate analysis. Bold indicates significant P-values (P<0.05).

Discussion

Identification of highly sensitive and specific biomarkers for prediction of outcome for patients with non-muscle invasive bladder tumours is of utmost clinical importance. Here we report the required step in clinical implementation of gene-expression signatures; the generation of a 12-gene signature for disease progression using qRT–PCR. We showed that the final 12-gene signature correlated highly significantly with outcome, and showed independent prognostic value when stratifying for stage, grade and treatment. Furthermore, by mutational analysis of genes frequently mutated in bladder cancer (FGFR3, PIK3CA and RAS family genes) combined with microarray-based analysis of an earlier reported 88-gene prognostic signature, we showed that intra-patient reproducibility and switch in tumour clonality may indicate that the 12-gene PCR signature analysis of several recurrent tumours is necessary in each single patient.

The analysis of intra-patient reproducibility of the microarray-based 88-gene molecular signature showed that in most cases (71%), identical classifications were obtained. In two of the five cases where divergent results were obtained, we also observed differences in the mutation status of FGFR3, indicating that the tumours may have arisen from different subclones. The monoclonal and oligoclonal nature of bladder tumours has been observed in several previous reports (Hartmann ; Zieger ). Sometimes the latest removed tumour was classified as high risk – sometimes it was the first removed tumour. This seems controversial; however, previous publications have indicated that the order of tumour resection may not be the order of formation of the tumour in the bladder. Very early tumours that grow very slowly may be resected later than offspring from this tumour that has a more rapid proliferation (van Tilborg ). Whether a switch in clonality is reflected at cystoscopic examination is at present unknown, but could be triggering a new molecular analysis. The analysis of intra-patient, and probably intra-tumour, variation of molecular signatures may be highly relevant in other cancer diseases with synchronous multifocal presentation like, for example, breast and prostate cancer.

Importantly, most of the samples used in this study to delineate the 12-gene PCR signature were analysed in the previous microarray-based validation study also, and consequently this study is to a large degree a technology-transfer study, and not a true validation of the molecular signature. Therefore, the relatively high HRs reported here for the 12-gene signature are based on the optimal cut-off values, and hence represent best case values, and the proposed 12-gene signature needs external validation before eventual clinical implementation. Furthermore, the technical reproducibility of the 12-gene signature needs to be further studied; however, because there is no need for normalisation genes, the assay and associated statistical calculations are easily handled.

The population studied here is neither consecutive nor included within a certain time interval, and consequently different biases may have been introduced; different follow-up guidelines and treatment regimens have been applied in different countries and in different time periods (e.g., a large fraction of high-risk patients have not received adjuvant treatment, and diagnosis of CIS from selected site biopsies has not been performed routinely). Furthermore, few women were included and several parameters regarding, for example, smoking history were not available. Therefore, prospective studies are needed for further validation of the findings. To address this, we are currently validating the 12-gene signature in an ongoing prospective study that will ultimately include 1200 patients with bladder cancer (UROMOL, FP7 EU project).

The 12-gene signature selected in this work includes six novel prognostic markers not included in previously reported signatures from our group. UBE2C and COL18A1 have been identified earlier in another meta-analysis of bladder cancer gene-expression data (Oncomine). The UBE2C and COL18A1 expression was described as being upregulated in progressing tumours and both genes were included in a 57-gene qRT–PCR signature (Wang ). SKAP2 has not been described in bladder cancer, but was earlier identified as having pro-invasion activities in human melanoma cells (Kabbarah ). NEK1 expression was identified in this study to be downregulated in progressing tumours compared with non-progressing tumours. Recently, this molecule has been shown to be involved in maintenance of chromosomal stability, and dysregulation of the molecule was associated with chromosomal instability (Chen ). COL4A3BP expression has been found increased in drug-resistant cell lines, and it has been suggested to be a target for chemotherapy-resistant cancers (Swanton ).

In conclusion, we have delineated a PCR-based 12-gene signature for disease progression in patients with non-muscle invasive bladder cancer. Based on molecular analysis of intra-patient variation, we suggest that this PCR-based test should be applied to tumours from all patient visits in order to increase test sensitivity.