Gene Expression Signature in Urine for Diagnosing and Assessing Aggressiveness of Bladder Urothelial Carcinoma

Cystoscopy and cytology are the standard methods used to detect and monitor bladder urothelial cell carcinoma (UCC). Cystoscopy is an invasive technique that has a high sensitivity (91%; ref. 1), except in cases of flat malignancies such as Tis. Cytology has the advantage of being noninvasive with a high specificity (90-96%; ref. 2), but it lacks sensitivity (11-76%; ref. 3), especially for low-grade disease.

Significant efforts have been made in developing noninvasive methods for detecting and predicting the biological behavior of UCC. In fact, the Food and Drug Administration has already accepted some of these tumor marker tests for use in routine patient care (4). Even though initial studies with these markers were promising, later reports often fail to confirm the results (5). It is likely that a panel of markers will overcome the limitations of these single markers.

The advent of microarray technology has had an important effect on the discovery of new molecular markers for diagnosing and predicting disease outcome in various cancers, including bladder cancer (6–12). Although gene expression has been extensively profiled in bladder tumor tissue samples, little is known about its correlation with gene expression profiles of urine samples. To this end, using microarrays, we identified a high number of differentially expressed genes between UCC and control samples as well as between low-grade (LG) and high-grade (HG) bladder tumors, in accordance with many other studies (6, 7, 13, 14). But more interestingly, we showed that these differentially expressed genes could also be identified in urine samples by real-time quantitative reverse-transcription PCR (qRT-PCR; ref. 12). Therefore, many of the differentially expressed genes detected could be candidates for the development of a noninvasive method of bladder cancer diagnosis and prediction of tumor aggressiveness.

The recently introduced approach of TaqMan Arrays (TA), based on the qRT-PCR platform, have already allowed the validation of large series of differentially expressed genes in bladder cancer (15). With the objective of achieving a noninvasive diagnostic and prognostic tool for UCC and following the hypothesis that a combination of genes could overcome the current single marker limitations, in this work we have analyzed the differentially expressed genes identified previously (12) in a large series of urine samples by means of TA methods.

The study is divided into two phases. In the first phase, a gene set expression signature has been obtained. During this phase, three consecutive steps were done wherein the number of genes analyzed progressively decreased while the number of samples evaluated was progressively increased (Fig. 1). In the second phase, we translated the gene set panel obtained to a routine clinical scenario by testing it in an independent series of voided urine samples.

A total of 803 urine samples were collected for the present study. Samples for phase I were collected between October 2003 and December 2004 in Fundació Puigvert, and samples for phase II were collected between January 2008 and July 2008 in Hospital Clínic of Barcelona. This study was approved by the Institutional Review Board of the two hospitals, and written informed consent to participate in the study was obtained from each patient and control before being enrolled.

After excluding those samples that did not fulfill the inclusion criteria (see Supplementary Data), the number of available samples for analysis was 576 (Table 1). In the first phase, 365 samples were analyzed: 244 urine samples from bladder washings (BW) from patients with histologically confirmed tumors (205 men, 39 women; mean age, 71.4 years) and 121 control voided urine samples (73 men, 48 women; mean age, 55.6 years). These samples were consecutively analyzed by qRT-PCR in TA: 39 were randomly chosen for analysis in TA that contained 384 genes (named 384-TA) selected from microarray data; 35 for analysis in TA that contained 96 genes (named 96-TA) selected from the former 384; and 291 for analysis in TA that contained a final group of 48 genes (named 48-TA; Table 2 and Supplementary Data). The grade and stage of the tumors were determined according to WHO criteria (16) and tumor-node-metastasis classification (17), respectively.

To compare gene expression values in BW and voided urine samples, 30 paired urine-BW samples from the same patient were collected from UCC patients. Of these, 3 paired samples were analyzed in 384-TA, 12 in 96-TA, and 15 in 48-TA. None of them was preamplified.

In the second phase, 211 voided urine samples were analyzed: 97 voided urine samples from UCC patients with histologically confirmed tumors (80 men, 17 women, mean age 72.6 years) and 114 control voided urine samples (90 men, 24 women, mean age 65.5 years; Table 1). These samples were analyzed in TAs that contained 48 genes (named 48-TA-urine).

BWs and voided urine samples were obtained and processed as described by Mengual et al. (18, 19). RNAs were extracted using TRIzol reagent (Invitrogen) according to the manufacturer's instructions. Total RNA was quantified by spectrophotometric analysis at 260 nm.

cDNA was synthesized from 100 ng to 1 μg of total RNA, depending on availability, using the High Capacity cDNA Archive Kit (Applied Biosystems, hereafter referred to as AB) following the manufacturer's instructions, except that the final volume of the reaction was 50 μL.

Multiplex preamplification of the 47 targets analyzed in the second phase (48-TA-urine) was achieved in cDNA from voided urine samples and was done as described by Mengual et al. (19), except that 12.5 μL of each cDNA sample were added to each reaction and that 14 preamplification cycles were done. 18S rRNA was not evaluated in 48-TA-urine because it cannot be included in the cDNA preamplification mix.

Target cDNAs were amplified in singleplex reactions using TA following the manufacturer's recommendations (AB). Briefly, 10 μL of cDNAs or 5 μL of preamplified cDNAs in the first and second part of the study, respectively, were mixed with 50 μL of 2× TaqMan Universal PCR MasterMix (AB) in a final volume of 100 μL. TAs were processed as described by Mengual et al. (19).

qRT-PCR data were processed with SDS 2.3 and RQ manager software packages (AB). For non-preamplified and preamplified cDNA samples, defined thresholds of 0.35 and 0.2, respectively, were used for all the genes to record the cycle thresholds (CT).

Relative expression level of target genes within a sample was expressed as ΔCT (ΔCT = CT_{mean of reference genes} − CT_{target gene}). Nonparametric Mann-Whitney test was done to examine statistically different expression patterns between groups. A P value of <0.05 was considered statistically significant. Heat maps were constructed using the heat map function under the R statistical environment (http://www.r-project.org/, http://sekhon.berkeley.edu/stats/html/heatmap.html). To cluster genes, a Pearson distance was calculated.

To correlate the 30 paired urine-BW samples, linear correlation coefficient of Pearson (R) was calculated.

The method to select the optimal and minimal combination of genes for diagnosis and aggressiveness assessment purposes from 48-TA has been a heuristic combination of median log 2 ratio, maximum independence among genes (using the Pearson correlation coefficient), and the visual analysis of frequency histograms. To combine the selected genes, we used a Naive Bayes classifier method (20). Briefly, this method consists in obtaining from every gene the sample distribution of expression frequencies in tumors and controls (which could be visualized, if data is categorized, as a histogram). Then, from this sample distribution, it is estimated the probability of the sample to be a tumor (T) if the expression of gene i (G_i) is equal to x or P(T|G_i = x), developing the Bayes theorem:

Thus, it is necessary to estimate, from the distribution of frequencies, P(G_i = x|T) and P(G_i = x|C) to obtain P(T|G_i = x).

Due to the continuous nature of data, it is necessary to establish some interval centered on the observed expression of G_i to estimate the required probabilities from the reference (or training) population. Initially, all samples from the reference population with expression values inside a small interval (0.05 around the sample observed expression value) are obtained. While the number of samples inside the interval (sum of tumors and controls) is lower than 15, the procedure is repeated recursively, increasing the interval by 0.05 on every step.

If independence is supposed between each gene, Bayes theorem could be applied recursively for each sample in all selected genes; thus, it is possible to use every a posteriori probability to be tumor (P(T|G_i = x) as a priori probability for the next gene [P(T|G_i = x) is assumed to be P(T) to calculate P(T|G_{i +1} = y) on the next Bayes iteration]. After the recursive application of Bayes theorem on all the selected genes, a final a posteriori probability is obtained for each sample. In our model, it is not required to establish a cutoff point to classify samples. If the a posteriori final probability is higher than the initial a priori probability [with an irrelevant initial value, in our case P(T) = 0.5 and P(C) = 0.5], the sample is classified as tumor. If it is lower, it is classified as control.

Because the samples of phase I were different from those of phase II, we do not have a test population for our final model. Consequently, we decided to use a leave-one-out cross-validation (LOOCV) technique to try to estimate the overfitting of the model, thus recalculating and comparing sensitivities and specificities for Naive Bayes algorithm (21).

To compare the performance of the Bayes classifier with that of another classification method, the same group of selected genes was analyzed using logistic regression. LOOCV technique was also used to estimate the overfitting of the logistic regression algorithm.

Normalized qRT-PCR data (ΔCT) were used to construct the heat map representations of cluster of the most differentially expressed genes (those with P < 0.05 and median log 2 ratio > ±1), grouped by sample type, in 384-TA and 96-TA. These heat maps showed different gene expression patterns in urine samples from control and from tumors as well as between LG and HG tumors. Heat maps with all the 45 finally selected genes for all cumulative samples (n = 365) showed 31 genes overexpressed and 14 genes underexpressed in tumor samples (Fig. 2). The differentially expressed genes (P < 0.05) selected to be analyzed by TA in each of the steps are listed in Supplementary Table S1.

According to the Bayes theorem, the final set of selected genes that provides the best classification between UCC and control samples was ANXA10, AHNAK2, CTSE, CRH, IGF2, KLF9, KRT20, MAGEA3, POSTN, PPP1R14D, SLC1A6, and TERT. Based on these genes, we obtained an overall sensitivity (SN) of 98% and specificity (SP) of 99% for the detection of bladder cancer using urine from BW. These values of SN and SP were maintained when grouping samples according to their pathologic stage (Table 3). Based on the same 12-gene combination, a classic logistic regression was calculated with similar results to the Bayes classification system (SN = 98% and SP = 94%). When the LOOCV method was used for Bayes system, a SN = 97% and a SP = 92% were obtained, whereas a SN = 96% and a SP = 93% were achieved for the logistic regression system, which indicates a general low overfitting effect.

Furthermore, we found that the addition of two genes (ASAM and MCM10) to the 12-gene set panel described for the diagnosis of UCC provided an overall SN = 79% and a SP = 92% to differentiate urine from patients with LG and HG tumors (Table 4). Based on this combination of genes, the logistic regression achieves a SN = 77% and a SP = 77%. When LOOCV was calculated for Bayes system, a SN = 68% and a SP = 77% were obtained, whereas for logistic regression the results were SN = 75% and SP = 75%. Although Naive Bayes presented a higher level of overfitting, the final performance of both techniques was similar.

The biological processes in which the 12+2 genes of the diagnostic panel are involved are summarized in Supplementary Data.

To compare the performance of BW and voided urine as samples for the diagnosis/prognosis of UCC, three paired BW-urine samples were initially run in 384-TA. As expected, the CT geometric mean of the three endogenous controls showed that urine RNA was more degraded than BW RNA. Therefore, the number of quantifiable genes (those genes with a CT value lower than 35) was lower in urine than in BW samples. However, high correlation coefficients (median R = 0.91) were obtained in these three samples when comparing gene expression values in BW and urine samples, indicating a high level of similarity between both samples (Supplementary Table S2). To confirm the findings found with 384-TA, 12 additional paired BW-urine samples were analyzed in 96-TA. The correlation between the measurements of quantified BW and urine target genes was comparable (median R = 0.87) with 384-TA. Finally, 15 additional paired BW-urine samples were analyzed in 48-TA, again obtaining comparable results (median R = 0.90).

At last, to determine the effect of analyzing BW instead of urine samples in the 12 diagnostic gene signature described here, we analyzed their performance in diagnosing bladder cancer in the 30 paired tumoral samples. Five of 30 cases were misclassified (one case was classified as control in both urine and BW and four cases were taken as controls when using voiding urine). These results raise the need of testing the gene signature in a new cohort of patients with voided urine. Clearly, the RNA obtained from urine samples is of lower quality than the RNA obtained from BW, suggesting that a previous preamplification of the sample can improve the consistency of the results.

Voided urine samples (n = 211) were used for the testing of the 12+2 gene set signature defined using BW samples. Because the effect of RNA degradation is much more apparent in this type of sample, a cDNA preamplification step was applied. All 48 genes from the third set of TA were tested again in voided urine samples. Figure 3 shows the heat map representations of all genes grouped by sample type and clustered by genes.

As expected, when testing the new preamplified voided urine samples from phase II with the 12-gene set model generated in the phase I using non-preamplified BW population, we found an important decrease in overall SN for the detection of bladder cancer in voided urine (SN = 70%) compared with BW samples while maintaining its specificity (SP = 96%; Table 3). We also found an overall decrement of SN and SP (SN = 42%; SP = 70%) in discriminating LG and HG tumors (Table 4).

To further refine the model in voided urine samples, we reanalyzed gene expressions of the 45 genes in urine samples (48-TA-urine) to include or exclude genes and to improve and adjust the gene expression signature. Interestingly, we found that the 12+2 gene set signature previously defined had the best performance in discriminating samples. Obviously, the estimation of probabilities for each gene was modified with the new series of data. Overall, SN and SP to differentiate urine from UCC patients and controls of the clinically applicable model were 89% and 95%, respectively. SN slightly decreased in the low-risk nonmuscle invasive bladder cancer (NMIBC) group (SN = 80%) but was higher in the high-risk NMIBC (SN = 93%) and muscle invasive bladder cancer (MIBC) groups (SN = 100%; Table 3). When a classic logistic regression was calculated for the same 12-gene combination, similar results to the Bayes classification system were obtained (SN = 86% and SP = 91%). After applying LOOCV for Bayes system, a SN = 81% and a SP = 91% were obtained, whereas for logistic regression the results of SN and SP were 84% and 89%, respectively. Again, there was more overfitting for the Bayes model, but the final performance of the model after LOOCV was very similar to logistic regression.

The overall SN and SP for discriminating urine from UCC patients with HG and LG tumors were 79% and 91%, respectively (Table 4). Based on this combination of genes, the logistic regression achieves a SN = 79% and a SP = 73%. When LOOCV was calculated for Bayes system, a SN = 68% and a SP = 68% were obtained, whereas for logistic regression the results were SN = 70% and SP = 61%. In this case, the Naive Bayes final prediction power seemed slightly superior than that of logistic regression.

Development of an accurate and noninvasive bladder tumor marker panel may be of great utility for screening, monitoring recurrence, and assessment of the true malignant potential of UCC. Thereby the number of surveillance cystoscopies done each year could be reduced and those patients with more aggressive tumors could be prioritized for treatment. Gene expression profiling has been used for the molecular classification of several cancers, including bladder cancer (6, 7, 12). However, to our knowledge, nobody has used voided urine samples to profile bladder cancer. In this study, gene expression profiles of urine are used to develop a noninvasive diagnostic and prognostic tool for bladder cancer. Our experimental approach was based on initially testing a few samples with a large number of previously selected genes (12) and progressively increasing the number of samples while reducing the number of genes. Finally, we translated the described gene set panel to a routine clinical scenario by testing it in an independent series of voided urine samples.

One of the major drawbacks for the clinical use of noninvasive RNA-based techniques to detect bladder cancer is the difficulty of obtaining sufficient quantities of high-quality RNA from urine samples. Our initial screenings showed that BW samples yielded higher amounts of and better-quality RNA than voided urine samples. Taking into account that when this study started, multiplex cDNA preamplification techniques were not available yet, we decided to use a minimally invasive obtained sample, the BW as a urine sample for tumor gene expression analysis. Because carrying out a BW in the case of controls could be considered an invasive maneuver, we used voided urine samples in these cases. Actually, gene expression comparison between 30 paired BW-urine tumoral samples using TA showed high correlation coefficients (median R = 0.90) on quantifiable genes in both samples, indicating a high degree of similarity between them, although differences are found. In addition, the number of quantifiable genes was lower in urine than in BW samples, suggesting that a cDNA preamplification technique could overcome this limitation.

The experimental design undertaken in this work to obtain the gene set signature was based on selecting the genes with the greatest differential expression in each set of TA. At the end of the selection steps, gene expression data for 45 genes were combined in all the samples (n = 244 tumors; n = 121 controls) to identify the gene set panel that best grouped UCC samples using the lowest number of genes. We found that a combination of 12 genes had the best performance in classifying tumor and control urine samples. We also found that the addition of gene expression information data from two additional genes to this panel provided information about tumor aggressiveness. The use of Bayes theorem to achieve the final gene set panel has worked slightly better than conventional logistic regression. Moreover, LOOCV has confirmed the consistency of the results.

Simultaneously with the course of this project, a multiplex cDNA preamplification technique and its applicability to TA technology (19) were developed. This has greatly facilitated our final purpose, which is to apply our gene set panel to voided urine samples. However, being aware that methodologic modifications could vary some of our results, we decided to analyze the complete 45-gene set again. The results showed that the 12+2 gene set panel described initially was the best in differentiating tumor-control and HG-LG tumors in voided urine samples. Although LOOCV indicates a certain degree of overfitting, all data obtained after cross-validation corroborate high SN and SP for the final model.

Urine cytology, the standard noninvasive method for detection and monitoring UCC patients, has a SP >90% in most studies (2). However, its SN ranges from 11% to 76% depending on tumor grade (3). This rather low SN is overcome with the use of cystoscopy, an invasive technique that shows an overall SN of 91% (1). Several noninvasive tumor markers have been developed with the aim of improving cytology results. Most of these Food and Drug Administration–approved markers have shown an increased SN with respect to cytology (58-79%), but none of them maintained its SP (66-87%; refs. 2, 5), UroVysion being the best performer (22). Taking into account the great heterogeneity of bladder tumors, we hypothesized that using a combination of different markers, we would achieve greater accuracy. Our test presents a SN and a SP in the high range (SN = 89%, SP = 95%), improving the results of currently available noninvasive UCC markers and tests. In addition, a qRT-PCR test has the advantage of being objective, easy to perform, and reliable.

Although high SN and SP levels are critical for any diagnostic/prognostic test, their applicability to the maximum number of samples is also essential. In this sense, only 10% of the voided urine samples are discarded in the second phase of this work because of low RNA quantity (<100 ng). However, we do not rule out the possibility that samples that yield less RNA quantity can be analyzed by TA after a cDNA preamplification process. Consequently, this percentage should be considered the maximum one, and we expect not to discard samples by this reason and to analyze virtually all samples in future studies.

Only two studies with a similar approach to ours are found in the literature. Holyoake et al. (23), using qRT-PCR technique, defined a combination of four individual markers to detect UCC in urine samples with a high SN for HG and MIBC samples (SN = 80-100%; SP = 85%), but with decreasing accuracy in LG and NMIBC samples (SN = 47-90%; SP = 85%). This is a phenomenon widely described for many UCC markers (3, 22). This trend is confirmed in our test, but achieving a SN of 80% in LG NMIBC, 93% in HG MIBC, and 100% in MIBC, while maintaining a SP of 95% in all three groups. It is noteworthy that this high SP has been achieved using not only healthy controls but also patients with nonneoplastic urologic diseases. By applying gene expression microarray analysis to cells recovered from BW, Rosser et al. (24) built a 14-gene molecular classifier that was able to classify patients with and without UCC with an overall accuracy of 76%. Surprisingly, none of the 319 differentially expressed gene probes found in this analysis is present in our final gene set expression signature. The use of different techniques, the limited number of patients analyzed, and the generation of the data from BW samples could account for this fact.

We also show in the present work that the gene expression signature of urine samples can provide information about disease aggressiveness at the time of initial diagnosis. The final model including two additional genes can differentiate LG from HG tumors with a SN of 79% and a SP of 91%. The assessment of the true malignant potential of a UCC at the time of diagnosis could change the current schedule of treatment of UCC. For instance, it could prioritize high-risk patients for an urgent cystoscopy or transurethral resection and could delay or prolong the interval between examinations in low-risk patients. Obviously, this greatly improves patient quality of life and outcome.

Limitations of this study include the lack of a final validation of the gene set expression panel using voided urine samples and a longer follow-up to predict patients outcome. In this sense, a multicentric randomized prospective study is currently being carried out, based on testing the gene set described, in a larger new population of urine samples from patients and controls. On the other hand, tumors can progress later than the time of follow-up in some patients. Thus, we cannot establish a gene set panel to predict outcome for patients with bladder cancer until a minimum follow-up of 5 years for all tumors has occurred.

In conclusion, this study describes that the gene expression signature present in urine is capable of identifying patients suffering from UCC. The utilization of such a sensitive and specific diagnostic/prognostic technique in clinical practice would improve detection and follow-up of UCC patients. However, introducing this noninvasive diagnosis/prognosis test for bladder cancer into clinical routine requires further validation in a prospective study using a larger number of voided urine samples.

No potential conflicts of interest were disclosed.

We thank all the staff from the Urology Departments and nurses from Fundació Puigvert and Hospital Clínic for collaborating in collecting urine samples. We are also grateful to pathology services from both hospitals for their support. We thank Helena Kruyer for the English correction of the manuscript, Dr. Dolors Colomer for assistance in the real-time PCR management, and Patricia Ruiz and Cristina Gazquez for technical support.

Grant Support: Laboratorios FINA BIOTECH. M. Ingelmo-Torres has a Research Support Staff Grant from Instituto Carlos III and M. Marín-Aguilera has a fellowship from IDIBAPS.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.