Abstract
Background: Bladder cancer is among the five most common malignancies worldwide, and due to high rates of recurrence, one of the most prevalent. Improvements in noninvasive urine-based assays to detect bladder cancer would benefit both patients and health care systems. In this study, the goal was to identify urothelial cell transcriptomic signatures associated with bladder cancer.
Methods: Gene expression profiling (Affymetrix U133 Plus 2.0 arrays) was applied to exfoliated urothelia obtained from a cohort of 92 subjects with known bladder disease status. Computational analyses identified candidate biomarkers of bladder cancer and an optimal predictive model was derived. Selected targets from the profiling analyses were monitored in an independent cohort of 81 subjects using quantitative real-time PCR (RT-PCR).
Results: Transcriptome profiling data analysis identified 52 genes associated with bladder cancer (P ≤ 0.001) and gene models that optimally predicted class label were derived. RT-PCR analysis of 48 selected targets in an independent cohort identified a 14-gene diagnostic signature that predicted the presence of bladder cancer with high accuracy.
Conclusions: Exfoliated urothelia sampling provides a robust analyte for the evaluation of patients with suspected bladder cancer. The refinement and validation of the multigene urothelial cell signatures identified in this preliminary study may lead to accurate, noninvasive assays for the detection of bladder cancer.
Impact: The development of an accurate, noninvasive bladder cancer detection assay would benefit both the patient and health care systems through better detection, monitoring, and control of disease. Cancer Epidemiol Biomarkers Prev; 21(12); 2149–58. ©2012 AACR.
This article is featured in Highlights of This Issue, p. 2123
Introduction
Bladder cancer is among the 5 most common malignancies worldwide (1). The current primary diagnostic approach to bladder cancer is cystoscopy coupled with voided urine cytology (VUC). Cystoscopy is an uncomfortable, invasive procedure associated with significant cost and possible infection and trauma. VUC remains the method of choice for the noninvasive detection of bladder cancer. However, while the assay has good specificity, VUC sensitivity is suboptimal, especially for low-grade and low-stage tumors (1, 2). A number of urine-based diagnostic protein markers have been developed commercially, but single biomarker assays lack adequate power to replace VUC and/or cystoscopy. This is not surprising given the redundancy of signaling pathways, the cross-talk between molecular networks, and the oligoclonality of tumors. Thus, it is postulated that the identification of alternative biomarkers for the detection of bladder cancer in urine may benefit from a more genome-wide analytic strategy.
In a previous pilot study of 46 cases (3), we showed the feasibility of gene expression profiling of exfoliated urothelia and developed an analytic approach to identify cancer-associated gene signatures. In this study, we have expanded the genome-wide expression profiling to 92 subjects (52 with confirmed bladder cancer). Parallel computational analyses were applied to identify bladder cancer-associated biomarkers and to derive multiplex diagnostic models. Selected candidate biomarkers were verified in an independent cohort of 81 subjects (44 bladder cancer cases) monitoring transcripts in urothelial cells obtained from naturally voided urine using quantitative real-time PCR (RT-PCR). An optimal diagnostic signature that comprised 14 genes achieved a high degree of accuracy. This discovery phase study shows the feasibility and potential of using urothelial cell gene expression signatures for the noninvasive detection of bladder cancer. Validation of the association of these models in larger and more diverse cohorts representing all urological disease states could lead to robust and accurate noninvasive tests for bladder cancer diagnosis.
Materials and Methods
Clinical sampling and processing
Under Institutional Review Board approval and informed consent, urine samples and associated clinical information were prospectively collected into a tissue bank. Two different clinical cohorts were analyzed in this study. The first group (cohort 1) consisted of 40 individuals with no previous history of urothelial cell carcinoma, gross hematuria, active urinary tract infection or urolithiasis (controls), and 52 individuals with newly diagnosed primary urothelia cell carcinoma. The second group (cohort 2) was comprised of 37 individuals with no previous history of urothelia cell carcinoma, gross hematuria, active urinary tract infection, or urolithiasis (controls), and 44 individuals with newly diagnosed primary urothelial cell carcinoma. All subjects were evaluated in the urology outpatient clinic, presenting with hematuria or voiding symptoms. All subjects underwent office cystoscopy, urinalysis, and VUC, and the majority also had axial imaging of the abdomen and pelvis. Furthermore, in our cancer group, histologic confirmation of urothelial cell carcinoma, including grade and stage, was noted from the excised tissue. In the discovery cohort, 50 mL of saline bladder barbotaged obtained during cystoscopy was collected into a sterile cup and stored at 4°C until it was processed (<3 hours). For the second cohort, 50 mL of midstream-voided urine was collected in a sterile cup and stored at 4°C until it was processed (<3 hours). Specimens with grossly visible blood were excluded. Pertinent information on clinical presentation and outcome were recorded. A summary of clinical data is given in Table 1. Each sample was assigned a unique identifying number before laboratory processing, as previously described (3). Laboratory personnel were blinded to final diagnosis.
. | Cohort 1 . | Cohort 2 . | ||
---|---|---|---|---|
. | Noncancer (%) . | Cancer (%) . | Noncancer (%) . | Cancer (%) . |
. | n = 40 . | n = 52 . | n = 37 . | n = 44 . |
Median age (range, y) | 69 (30–90) | 68 (36–90) | 69 (19–79) | 67 (47–90) |
Male/female ratio | 26:14 | 42:10 | 21:16 | 36:8 |
Race | ||||
White | 30 (75) | 48 (92) | 32 (86) | 40 (91) |
African American | 7 (18) | 3 (6) | 4 (11) | 3 (7) |
Other | 3 (7) | 1 (1) | 1 (3) | 1 (2) |
Tobacco history | 22 (55) | 46 (88) | 16 (43) | 33 (75) |
Suspicious/positive cytology | 0 (0) | 16 (31) | 0 (0) | 12 (27) |
Clinical stage | ||||
Tis | n/a | 3 (6) | n/a | 0 (0) |
Ta | n/a | 11 (21) | n/a | 9 (20) |
T1 | n/a | 16 (31) | n/a | 8 (18) |
T2 | n/a | 22 (42) | n/a | 24 (55) |
T3 | n/a | 0 (0) | n/a | 3 (7) |
Grade | ||||
1/2 | n/a | 9 (16) | n/a | 7 (16) |
3 | n/a | 43 (84) | n/a | 37 (84) |
. | Cohort 1 . | Cohort 2 . | ||
---|---|---|---|---|
. | Noncancer (%) . | Cancer (%) . | Noncancer (%) . | Cancer (%) . |
. | n = 40 . | n = 52 . | n = 37 . | n = 44 . |
Median age (range, y) | 69 (30–90) | 68 (36–90) | 69 (19–79) | 67 (47–90) |
Male/female ratio | 26:14 | 42:10 | 21:16 | 36:8 |
Race | ||||
White | 30 (75) | 48 (92) | 32 (86) | 40 (91) |
African American | 7 (18) | 3 (6) | 4 (11) | 3 (7) |
Other | 3 (7) | 1 (1) | 1 (3) | 1 (2) |
Tobacco history | 22 (55) | 46 (88) | 16 (43) | 33 (75) |
Suspicious/positive cytology | 0 (0) | 16 (31) | 0 (0) | 12 (27) |
Clinical stage | ||||
Tis | n/a | 3 (6) | n/a | 0 (0) |
Ta | n/a | 11 (21) | n/a | 9 (20) |
T1 | n/a | 16 (31) | n/a | 8 (18) |
T2 | n/a | 22 (42) | n/a | 24 (55) |
T3 | n/a | 0 (0) | n/a | 3 (7) |
Grade | ||||
1/2 | n/a | 9 (16) | n/a | 7 (16) |
3 | n/a | 43 (84) | n/a | 37 (84) |
Gene expression profiling
Gene expression profiling was conducted on Affymetrix Human Genome arrays (U133 Plus 2.0) according to standard protocols (Affymetrix). All microarray data obtained in the course of this study have been deposited in the NCBI's Gene Expression Omnibus (4) and are accessible through the GEO Series accession number GSE31189. Conventional comparative analysis of gene expression data for association of genes with bladder cancer was conducted as previously described (3, 5, 6). Briefly, the 2-sample Welch t statistics that allow unequal variances were used to identify genes that were differentially expressed between normal and tumor samples. The P value was used to assess the statistical significance for each gene. To correct for multitesting errors, the family-wise error rate (FWER) was used following a permutation-based bootstrap step-down minP procedure.
Derivation of a diagnostic gene signature from microarray data
To derive optimal multigene diagnostic signatures from the microarray data, we applied the LoGo feature selection algorithm (7, 8). To avoid possible overfitting of computational models to training data, we used the leave-one-out cross validation (LOOCV) method to estimate classifier parameters and prediction performance (9). A receiver operating characteristic (ROC) curve (10), obtained by varying a decision threshold, was used to provide a direct view of how a prediction model functioned at different sensitivity and specificity levels, and significance between groups was evaluated using t tests.
Quantitative RT-PCR analysis
Purified RNA samples were evaluated quantitatively and qualitatively using an Agilent Bioanalyzer 2000. Complementary DNA was synthesized from 20 to 500 ng of total RNA, using the High Capacity cDNA Reverse Transcriptase Kit (Applied Biosystems) according the manufacturer's instructions, with random primers in a total reaction volume of 20 μL.
Selection of endogenous reference controls.
An aliquot of each sample cDNA was used in a multiplex RT-PCR preamplification reaction of 15 endogenous reference targets: GAPDH, ACTB, B2M, GUSB, HMBS, HPRT1, IPO8, PGK1, POLR2A, PPIA, RPLP0, TBP, TFRC, UBC, and YWHAZ. The 15 TaqMan Gene Expression Assays were pooled together at 0.2× final concentration. Subsequently, 12.5 μL of the pooled assay mix (0.2×) was combined with 4 μL of each cDNA sample and 25 μL of the TaqMan PreAmp Master Mix (2×) in a final volume of 50 μL. Thermal cycling conditions were as follows: initial hold at 95°C during 10 minutes and 10 preamplification cycles of 15 seconds at 95°C and 4 minutes at 60°C. The preamplification products were diluted 1:5 with Tris-EDTA buffer before singleplex reaction amplification using the TaqMan Human Endogenous Control Array (PN 4367563; Applied Biosystems), a 384-well microfluidic card containing the 15 genes listed above plus 18S RNA. The reactions were carried out on a 7900HT Fast Real-Time PCR System (Applied Biosystems). Genes with the least variable expression across all 81 samples (UBC, PPIA, PGK1, and GAPDH) were identified using GeNorm software (Integromics) and subsequently selected for the 48-target custom TaqMan Low Density Array (TLDA). The TLDA format is a 384-well system that uses standard TaqMan assays and enables automated loading and high-throughput analyses (11).
Custom array preamplification and amplification reactions.
TLDAs were constructed by Applied Biosystems using predesigned assays whose probe would span exons. Targets included were: UBC, PPIA, PGK1, and GAPDH (4 endogenous controls); BIRC5, TERT, KRT20, CLU, PLAU, CALR, ANG, CA9, SCG3, ATF3, TLR2, AGT, DMBT1, ERBB2, CTNNA1, ATM, FBXO9, CCNE2, NUSAP1, SNAI2, PLOD2, MMP12, IL1RAP, ITGB5, DSC2, APOE, TMEM45A, SYNGR1, MMP10, IL8, VEGFA, CPA1, CCL18, CRH, MCOLN1 SERPINE1 MMP1, MMP9, FGF22, MXRA8, NRG3, SEMA3D, PTX3, and RAB1A (44 biomarker targets). A multiplex RT-PCR preamplification reaction was carried out using the pooled 48 TaqMan Gene Expression Assays. Assay reagents at 0.2× final concentration were combined with 7.5 μL of each cDNA sample and 15 μL of the TaqMan PreAmp Master Mix (2×) in a final volume of 30 μL. Thermal cycling conditions were as follows: initial hold at 95°C during 10 minutes; 14 preamplification cycles of 15 seconds at 95°C, 4 minutes at 60°C, and a final hold at 99.9°C for 10 minutes. Ten microliters of undiluted preamplification products was used in the subsequent singleplex amplification reactions, combined with 50 μL of 2× TaqMan Universal PCR MasterMix (Applied Biosystems) in a final volume of 100 μL, according to the manufacturer's instructions. One sample of Human Universal Reference Total cDNA (Clontech) was included as a calibrator in each microfluidic card. The reactions were run in a 7900HT Fast Real-Time PCR System (Applied Biosystems).
Statistical analysis
RT-PCR amplification results were processed with RQ manager (Applied Biosystems) and StatMiner (Integromics) software packages. The baseline correction was manually checked for each target and the Ct threshold was set to 0.2 for every target across all plates. One sample that amplified less than 10% of the 48 targets was removed from the analysis, and 1 target (PTX3) was removed from all subsequent analyses because amplification curve analysis revealed amplification of nonspecific products. The RT-PCR data was analyzed using the comparative Ct method (12). ΔΔCt values were calculated using a geometric average of the 4 endogenous reference targets (UBC, PPIA, PGK1, and GAPDH) as normalizer and Human Universal Reference Total cDNA (Clontech) as the calibrator. Genes deemed to be differentially expressed between bladder cancer and noncancer samples were determined by t test comparison (P < 0.01).
Logistic regression was used to establish a prediction model from the quantitative PCR data, that is, to predict the actual status of a given sample as cancer or control (see Supplementary Document SA). A ROC curve was then plotted to visualize how a prediction model carried out at all sensitivity and specificity levels, and the area under ROC curves (AUC) was reported following the standards for the reporting of diagnostic accuracy studies criteria for biomarker study reporting (13). The P value was computed as the occurrence frequencies of the iterations, where the resulting AUCs outperformed that obtained using the original class labels. A P < 0.01 was considered to be statistically significant. Logistic regression analysis was conducted with bladder cancer status (yes vs. no) as the response variable and biomarker values, gender, and ethnicity. Neither gender nor ethnicity was identified as a confounder. We also conducted a permutation test (1000-fold) to estimate the P value of predictive performance. Further details of model derivation are available in Supplementary document SA. Statistical analyses were conducted by SPSS 13.0. (IBM) and by MedCalc version 8.0 (MedCalc Software).
Results
The overall scheme for identifying diagnostic molecular signatures for bladder cancer is depicted in Fig. 1. Urothelial cell samples were isolated from 92 subjects for molecular profile analysis (discovery cohort). Of these subjects, 52 had biopsy confirmed bladder cancer and 40 had no evidence of cancer. Patient cohort characteristics are summarized in Table 1.
Gene expression profiling
RNA was subjected to 2-cycle amplification strategy before hybridization to Affymetrix U133 Plus 2.0 arrays. Comparative group (cancer vs. noncancer) data analysis identified a 447-gene set that had expression patterns associated with disease status (P ≤ 0.01), with 52 genes having a P value of 0.001 or less (see Supplementary Table S1). From a list of genes ranked only by P value, it is not clear which are the most relevant to the classification task at hand, in this case, stratifying cancer versus noncancer cases. To identify the most accurate diagnostic signatures from the microarray data, we applied a feature selection algorithm that we have previously derived and applied to the derivation of optimal disease classifiers in bladder, breast, and prostate cancer (3, 7–9). The algorithm conducts multivariate data analyses on high-dimensional data without making any assumptions about the underlying data distribution. This approach identified a 43-gene model (Supplementary Table S2) that functioned best in predicting class label, achieving an area under the ROC curve of 0.82, during leave-one-out cross validation. Overlap between the modeling and comparative data analyses informed the selection of genes for investigation with quantitative PCR in the second cohort.
Verification of candidate biomarker transcripts in an independent cohort
To validate and refine the bladder cancer diagnostic molecular signature, we tested a selected panel of candidate biomarkers from the analyses described above in naturally voided urine samples from an independent cohort (cohort 2; Table 1) composed of 44 subjects with biopsy-proven bladder cancer and 37 subjects with no evidence of cancer (Table 1). Target transcripts were measured in RNA samples from urothelia isolated from naturally micturated urine using quantitative RT-PCR. Microfluidic cards of 384-well format were custom designed to include 44 candidate biomarker targets, plus 4 selected endogenous controls selected by screening the level of 15 commonly used endogenous controls in the full cohort of samples (selection and design described in Methods and in Supplementary document SB). Biomarker targets were selected primarily from the top statistical ranking and molecular signature models described above, but we also included several putative biomarkers (TERT, KRT20, CLU, PLAU, CALR, CA9, and ANG) that have been shown to be associated with bladder cancer in our own previous studies or the literature (14, 15). When other selection criteria were equal, we selected genes that encode integral membrane proteins or secreted proteins because these classes hold particular potential for development as urine-based biomarker assays. The 44 selected targets are listed in Table 2. Differential expression values were calculated by normalization using the reference targets (UBC, PPIA, PGK1, and GAPDH) and Human Universal Reference Total cDNA (Clontech) as the calibrator on each plate. Fifteen targets were revealed to be differentially expressed between bladder cancer and noncancer samples (P < 0.05) using a t test.
Gene symbol . | Gene name . | Gene ID . | Mean ΔΔCt N . | SD ΔΔCt N . | Mean ΔΔCt bladder cancer . | SD ΔΔCt bladder cancer . | t test P value . | FCa bladder cancer/N . |
---|---|---|---|---|---|---|---|---|
CA9 | carbonic anhydrase IX | 768 | 9.11 | 5.25 | 2.07 | 4.53 | 1.E-08 | 131.63 |
MMP12 | matrix metallopeptidase 12 | 4321 | 1.20 | 5.05 | −4.27 | 4.08 | 9.E-07 | 44.39 |
CCL18 | chemokine (c-c motif) ligand 18 | 6362 | 5.76 | 4.96 | 1.30 | 3.00 | 5.E-06 | 22.03 |
MMP10 | matrix metallopeptidase 10 | 4319 | 4.69 | 5.16 | −0.30 | 5.05 | 4.E-05 | 31.75 |
TMEM45A | transmembrane protein 45A | 55076 | 6.62 | 5.68 | 2.13 | 3.93 | 9.E-05 | 22.43 |
ANG | angiogenin, ribonuclease, RNase A family, 5 | 283 | 9.64 | 4.31 | 13.02 | 3.96 | 0.001 | 0.1 |
SNAI2 | snail homolog 2 (Drosophila) | 6591 | 5.85 | 5.08 | 2.75 | 3.07 | 0.001 | 8.56 |
MMP1 | matrix metallopeptidase 1 | 4312 | 0.89 | 4.84 | −2.78 | 4.91 | 0.001 | 12.77 |
SERPINE1 | serpin peptidase inhibitor, clade E, member 1 | 5054 | 4.14 | 3.87 | 1.38 | 4.01 | 0.003 | 6.74 |
MXRA8 | matrix–remodeling-associated protein 8 | 54587 | 9.59 | 4.08 | 6.96 | 3.68 | 0.004 | 6.19 |
MMP9 | matrix metallopeptidase 9 | 4318 | −0.76 | 3.91 | −3.17 | 3.26 | 0.004 | 5.32 |
CCNE2 | cyclin E2 | 9134 | 3.52 | 4.45 | 1.06 | 4.14 | 0.013 | 5.5 |
BIRC5 | baculoviral IAP repeat containing 5 | 332 | 7.57 | 4.72 | 5.04 | 4.59 | 0.018 | 5.79 |
SEMA3D | semaphorin-3D | 223117 | 7.43 | 4.18 | 5.32 | 4.12 | 0.027 | 4.3 |
PLAU | plasminogen activator, urokinase | 5328 | 1.22 | 3.11 | −0.22 | 3.32 | 0.052 | 2.71 |
SYNGR1 | synaptogyrin 1 | 9145 | 0.44 | 3.65 | 1.97 | 3.77 | 0.073 | 0.35 |
SCG3 | secretogranin III | 29106 | 10.63 | 4.76 | 8.56 | 5.30 | 0.074 | 4.19 |
NUSAP1 | nucleolar and spindle associated protein 1 | 51203 | 1.74 | 3.97 | 0.31 | 3.53 | 0.096 | 2.68 |
CRH | corticotropin releasing hormone | 1392 | 9.12 | 4.44 | 6.97 | 6.84 | 0.109 | 4.43 |
TERT | telomerase reverse transcriptase | 7015 | 7.07 | 4.94 | 4.96 | 6.37 | 0.110 | 4.3 |
DMBT1 | deleted in malignant brain tumors 1 | 1755 | 2.91 | 5.76 | 5.01 | 6.37 | 0.131 | 0.23 |
IL1RAP | interleukin 1 receptor accessory protein | 3556 | -0.53 | 3.32 | −1.72 | 3.57 | 0.132 | 2.28 |
NRG3 | neuregulin 3 | 10718 | 8.26 | 5.15 | 9.74 | 4.95 | 0.195 | 0.36 |
RAB1A | RAB1A, member RAS oncogene family | 5861 | 1.58 | 1.64 | 2.24 | 2.74 | 0.211 | 0.63 |
AGT | angiotensinogen | 183 | 8.16 | 4.44 | 9.36 | 4.75 | 0.254 | 0.44 |
PLOD2 | procollagen-lysine, 2-oxoglutarate 5-dioxygenase 2 | 5352 | 4.51 | 4.25 | 3.41 | 4.27 | 0.254 | 2.15 |
MCOLN1 | mucolipin 1 | 57192 | 0.76 | 3.69 | −0.22 | 3.92 | 0.260 | 1.97 |
FBXO9 | f-box protein 9 | 26268 | 1.09 | 2.89 | 1.90 | 3.69 | 0.287 | 0.57 |
IL8 | interleukin 8 | 3576 | −6.26 | 4.03 | −7.55 | 6.58 | 0.308 | 2.44 |
CTNNA1 | catenin (cadherin-associated protein), alpha 1, 102kDa | 1495 | 0.53 | 2.27 | 1.17 | 3.24 | 0.319 | 0.64 |
CLU | Clusterin | 1191 | 5.51 | 2.99 | 6.22 | 3.59 | 0.347 | 0.61 |
ATM | ataxia telangiectasia mutated | 472 | 0.52 | 4.11 | −0.18 | 3.36 | 0.406 | 1.63 |
VEGFA | vascular endothelial growth factor A | 7422 | −0.50 | 2.48 | −0.95 | 4.45 | 0.590 | 1.37 |
CALR | Calreticulin | 811 | 1.33 | 2.78 | 0.91 | 4.12 | 0.599 | 1.34 |
DSC2 | desmocollin 2 | 1824 | 0.20 | 2.72 | 0.60 | 4.03 | 0.613 | 0.76 |
FGF22 | fibroblast growth factor 22 | 27006 | 2.36 | 3.68 | 1.96 | 3.51 | 0.623 | 1.32 |
CPA1 | carboxypeptidase A1 (pancreatic) | 1357 | 12.97 | 4.34 | 13.48 | 5.43 | 0.650 | 0.7 |
APOE | apolipoprotein E | 348 | 5.46 | 2.81 | 5.19 | 2.99 | 0.674 | 1.21 |
ATF3 | activating transcription factor 3 | 467 | −0.34 | 2.58 | −0.09 | 3.38 | 0.711 | 0.84 |
ERBB2 | v-erb-b2 erythroblastic leukemia viral oncogene homolog 2 | 2064 | −0.64 | 4.00 | −0.89 | 3.71 | 0.774 | 1.19 |
KRT20 | keratin 20 | 54474 | −10.19 | 4.46 | −10.45 | 5.53 | 0.816 | 1.2 |
TLR2 | toll-like receptor 2 | 7097 | −1.07 | 2.73 | −1.24 | 4.50 | 0.838 | 1.13 |
ITGB5 | integrin, beta 5 | 3693 | 2.70 | 3.39 | 2.72 | 3.24 | 0.977 | 0.99 |
UBC | ubiquitin C | 7316 | Endog. Control | NA | ||||
PPIA | peptidylprolyl isomerase A | 5478 | Endog. Control | NA | ||||
PGK1 | phosphoglycerate kinase 1 | 5230 | Endog. Control | NA | ||||
GAPDH | glyceraldehyde-3-phosphate dehydrogenase | 2597 | Endog. Control | NA |
Gene symbol . | Gene name . | Gene ID . | Mean ΔΔCt N . | SD ΔΔCt N . | Mean ΔΔCt bladder cancer . | SD ΔΔCt bladder cancer . | t test P value . | FCa bladder cancer/N . |
---|---|---|---|---|---|---|---|---|
CA9 | carbonic anhydrase IX | 768 | 9.11 | 5.25 | 2.07 | 4.53 | 1.E-08 | 131.63 |
MMP12 | matrix metallopeptidase 12 | 4321 | 1.20 | 5.05 | −4.27 | 4.08 | 9.E-07 | 44.39 |
CCL18 | chemokine (c-c motif) ligand 18 | 6362 | 5.76 | 4.96 | 1.30 | 3.00 | 5.E-06 | 22.03 |
MMP10 | matrix metallopeptidase 10 | 4319 | 4.69 | 5.16 | −0.30 | 5.05 | 4.E-05 | 31.75 |
TMEM45A | transmembrane protein 45A | 55076 | 6.62 | 5.68 | 2.13 | 3.93 | 9.E-05 | 22.43 |
ANG | angiogenin, ribonuclease, RNase A family, 5 | 283 | 9.64 | 4.31 | 13.02 | 3.96 | 0.001 | 0.1 |
SNAI2 | snail homolog 2 (Drosophila) | 6591 | 5.85 | 5.08 | 2.75 | 3.07 | 0.001 | 8.56 |
MMP1 | matrix metallopeptidase 1 | 4312 | 0.89 | 4.84 | −2.78 | 4.91 | 0.001 | 12.77 |
SERPINE1 | serpin peptidase inhibitor, clade E, member 1 | 5054 | 4.14 | 3.87 | 1.38 | 4.01 | 0.003 | 6.74 |
MXRA8 | matrix–remodeling-associated protein 8 | 54587 | 9.59 | 4.08 | 6.96 | 3.68 | 0.004 | 6.19 |
MMP9 | matrix metallopeptidase 9 | 4318 | −0.76 | 3.91 | −3.17 | 3.26 | 0.004 | 5.32 |
CCNE2 | cyclin E2 | 9134 | 3.52 | 4.45 | 1.06 | 4.14 | 0.013 | 5.5 |
BIRC5 | baculoviral IAP repeat containing 5 | 332 | 7.57 | 4.72 | 5.04 | 4.59 | 0.018 | 5.79 |
SEMA3D | semaphorin-3D | 223117 | 7.43 | 4.18 | 5.32 | 4.12 | 0.027 | 4.3 |
PLAU | plasminogen activator, urokinase | 5328 | 1.22 | 3.11 | −0.22 | 3.32 | 0.052 | 2.71 |
SYNGR1 | synaptogyrin 1 | 9145 | 0.44 | 3.65 | 1.97 | 3.77 | 0.073 | 0.35 |
SCG3 | secretogranin III | 29106 | 10.63 | 4.76 | 8.56 | 5.30 | 0.074 | 4.19 |
NUSAP1 | nucleolar and spindle associated protein 1 | 51203 | 1.74 | 3.97 | 0.31 | 3.53 | 0.096 | 2.68 |
CRH | corticotropin releasing hormone | 1392 | 9.12 | 4.44 | 6.97 | 6.84 | 0.109 | 4.43 |
TERT | telomerase reverse transcriptase | 7015 | 7.07 | 4.94 | 4.96 | 6.37 | 0.110 | 4.3 |
DMBT1 | deleted in malignant brain tumors 1 | 1755 | 2.91 | 5.76 | 5.01 | 6.37 | 0.131 | 0.23 |
IL1RAP | interleukin 1 receptor accessory protein | 3556 | -0.53 | 3.32 | −1.72 | 3.57 | 0.132 | 2.28 |
NRG3 | neuregulin 3 | 10718 | 8.26 | 5.15 | 9.74 | 4.95 | 0.195 | 0.36 |
RAB1A | RAB1A, member RAS oncogene family | 5861 | 1.58 | 1.64 | 2.24 | 2.74 | 0.211 | 0.63 |
AGT | angiotensinogen | 183 | 8.16 | 4.44 | 9.36 | 4.75 | 0.254 | 0.44 |
PLOD2 | procollagen-lysine, 2-oxoglutarate 5-dioxygenase 2 | 5352 | 4.51 | 4.25 | 3.41 | 4.27 | 0.254 | 2.15 |
MCOLN1 | mucolipin 1 | 57192 | 0.76 | 3.69 | −0.22 | 3.92 | 0.260 | 1.97 |
FBXO9 | f-box protein 9 | 26268 | 1.09 | 2.89 | 1.90 | 3.69 | 0.287 | 0.57 |
IL8 | interleukin 8 | 3576 | −6.26 | 4.03 | −7.55 | 6.58 | 0.308 | 2.44 |
CTNNA1 | catenin (cadherin-associated protein), alpha 1, 102kDa | 1495 | 0.53 | 2.27 | 1.17 | 3.24 | 0.319 | 0.64 |
CLU | Clusterin | 1191 | 5.51 | 2.99 | 6.22 | 3.59 | 0.347 | 0.61 |
ATM | ataxia telangiectasia mutated | 472 | 0.52 | 4.11 | −0.18 | 3.36 | 0.406 | 1.63 |
VEGFA | vascular endothelial growth factor A | 7422 | −0.50 | 2.48 | −0.95 | 4.45 | 0.590 | 1.37 |
CALR | Calreticulin | 811 | 1.33 | 2.78 | 0.91 | 4.12 | 0.599 | 1.34 |
DSC2 | desmocollin 2 | 1824 | 0.20 | 2.72 | 0.60 | 4.03 | 0.613 | 0.76 |
FGF22 | fibroblast growth factor 22 | 27006 | 2.36 | 3.68 | 1.96 | 3.51 | 0.623 | 1.32 |
CPA1 | carboxypeptidase A1 (pancreatic) | 1357 | 12.97 | 4.34 | 13.48 | 5.43 | 0.650 | 0.7 |
APOE | apolipoprotein E | 348 | 5.46 | 2.81 | 5.19 | 2.99 | 0.674 | 1.21 |
ATF3 | activating transcription factor 3 | 467 | −0.34 | 2.58 | −0.09 | 3.38 | 0.711 | 0.84 |
ERBB2 | v-erb-b2 erythroblastic leukemia viral oncogene homolog 2 | 2064 | −0.64 | 4.00 | −0.89 | 3.71 | 0.774 | 1.19 |
KRT20 | keratin 20 | 54474 | −10.19 | 4.46 | −10.45 | 5.53 | 0.816 | 1.2 |
TLR2 | toll-like receptor 2 | 7097 | −1.07 | 2.73 | −1.24 | 4.50 | 0.838 | 1.13 |
ITGB5 | integrin, beta 5 | 3693 | 2.70 | 3.39 | 2.72 | 3.24 | 0.977 | 0.99 |
UBC | ubiquitin C | 7316 | Endog. Control | NA | ||||
PPIA | peptidylprolyl isomerase A | 5478 | Endog. Control | NA | ||||
PGK1 | phosphoglycerate kinase 1 | 5230 | Endog. Control | NA | ||||
GAPDH | glyceraldehyde-3-phosphate dehydrogenase | 2597 | Endog. Control | NA |
NOTE: TaqMan Assay IDs for the 44 target genes and 4 endogenous controls are in Supplementary Table S3.
aFold change = |${\rm 2}^{- \Delta \Delta C_{\rm t}} {\rm T/2}^{- \Delta \Delta C_{\rm t}} {\rm N}$| (See Materials and Methods Section).
N, noncancer control
Molecular signature derivation
To identify the molecular signature that could best predict the status of a given sample as cancer or control, we also used L1-regularized logistic regression to derive a prediction model from the 44 targets. ROC curves were generated to compare the performance of all possible combinatorial models. The optimal signature derived from the PCR analysis was composed of 14 genes (Table 3). We observe that some genes (e.g., DMBT1 and ERBB2) have relatively high P values and, therefore, may not be valuable when evaluated individually; however, these targets do provide critical information when combined with a panel of biomarkers. A scatter plot of the 14-gene signature performance showed a significant (P < 0.0001) spread of values between cancer and noncancer cases (Fig. 2A), and the ROC curve depicts performance across all sensitivity and specificity values (Fig. 2B). The ROC curve reveals that, at a sensitivity of 90%, the 14-gene signature achieved a specificity of 100%. The AUC was 0.98. By plotting the variation of signature P values against AUC values, it was revealed that the 9 top-ranked genes contributed the most to the 14-gene signature performance (Fig. 2C).
Gene symbol . | Gene name . | Fold change bladder cancer/N . | 14-gene signature average weighta . | 9-gene signature average weighta . | Localization . |
---|---|---|---|---|---|
CA9 | carbonic anhydrase IX | 131.63 | 0.18 | 0.480 | Cell membrane/nucleus |
MMP12 | matrix metallopeptidase 12 | 44.39 | 0.075 | Secreted | |
TMEM45A | transmembrane protein 45A | 22.43 | 0.06 | 0.032 | Membrane |
CCL18 | chemokine (C-C motif) ligand 18 | 22.03 | 0.16 | 0.400 | Secreted |
MXRA8 | matrix-remodelling associated 8 | 6.19 | 0.02 | Membrane | |
MMP9 | matrix metallopeptidase 9 | 5.32 | 0.05 | 0.009 | Secreted |
CRH | corticotropin releasing hormone | 4.43 | 0.001 | Secreted | |
SEMA3D | semaphorin 3D | 4.30 | 0.05 | 0.002 | Secreted |
VEGFA | vascular endothelial growth factor A | 1.37 | −0.02 | Secreted | |
ERBB2 | v-erb-b2 erythroblastic leukemia viral oncogene homolog 2 | 1.19 | 0.08 | 0.001 | Membrane/cytoplasm nucleus |
DSC2 | desmocollin 2 | 0.76 | −0.04 | 0.001 | Cell membrane |
RAB1A | RAB1A, member RAS oncogene family | 0.63 | −0.05 | Endoplasmic reticulum | |
AGT | angiotensinogen | 0.44 | −0.01 | Secreted | |
SYNGR1 | synaptogyrin 1 | 0.35 | −0.06 | Membrane | |
DMBT1 | deleted in malignant brain tumors 1 | 0.23 | −0.09 | Secreted | |
ANG | angiogenin, ribonuclease, RNase A family, 5 | 0.10 | −0.08 | Secreted |
Gene symbol . | Gene name . | Fold change bladder cancer/N . | 14-gene signature average weighta . | 9-gene signature average weighta . | Localization . |
---|---|---|---|---|---|
CA9 | carbonic anhydrase IX | 131.63 | 0.18 | 0.480 | Cell membrane/nucleus |
MMP12 | matrix metallopeptidase 12 | 44.39 | 0.075 | Secreted | |
TMEM45A | transmembrane protein 45A | 22.43 | 0.06 | 0.032 | Membrane |
CCL18 | chemokine (C-C motif) ligand 18 | 22.03 | 0.16 | 0.400 | Secreted |
MXRA8 | matrix-remodelling associated 8 | 6.19 | 0.02 | Membrane | |
MMP9 | matrix metallopeptidase 9 | 5.32 | 0.05 | 0.009 | Secreted |
CRH | corticotropin releasing hormone | 4.43 | 0.001 | Secreted | |
SEMA3D | semaphorin 3D | 4.30 | 0.05 | 0.002 | Secreted |
VEGFA | vascular endothelial growth factor A | 1.37 | −0.02 | Secreted | |
ERBB2 | v-erb-b2 erythroblastic leukemia viral oncogene homolog 2 | 1.19 | 0.08 | 0.001 | Membrane/cytoplasm nucleus |
DSC2 | desmocollin 2 | 0.76 | −0.04 | 0.001 | Cell membrane |
RAB1A | RAB1A, member RAS oncogene family | 0.63 | −0.05 | Endoplasmic reticulum | |
AGT | angiotensinogen | 0.44 | −0.01 | Secreted | |
SYNGR1 | synaptogyrin 1 | 0.35 | −0.06 | Membrane | |
DMBT1 | deleted in malignant brain tumors 1 | 0.23 | −0.09 | Secreted | |
ANG | angiogenin, ribonuclease, RNase A family, 5 | 0.10 | −0.08 | Secreted |
Abbreviation: N, noncancer control
aGenes sorted in descending order of fold change. Negative prefix denotes a gene downregulated in cancer cases.
The optimal 14-gene model was composed of 7 gene transcripts that were upregulated in urines from bladder cancer subjects and 7 that were downregulated. Considering the development of potential biomarker assays, in which upregulated genes are more applicable for accurate detection, we applied a constrained L1-regularized logistic regression model (Supplementary document SA) to derive an optimal diagnostic signature composed of only upregulated genes (Supplementary document SA). Considering the limits on target selection, the 9-gene signature (Table 3) worked very well, achieving an AUC of 0.93 (Fig. 3A and B). From the ROC curve (Fig. 3), we can see that, at a sensitivity of 80%, the 9-gene signature achieved 98% specificity. Plotting the P value variation and AUC values revealed that the top 2 ranked genes (CA9 and CCL18) contributed most to the 9-gene signature (Fig. 3C).
Discussion
A number of molecular tests have been developed to detect potential tumor-associated biomarkers that are released into the urine from malignant cells. Some of these tests include bladder tumor antigen (BTA; ref. 16), nuclear matrix protein 22 (NMP-22; ref. 17), and FISH analysis for chromosomal abnormalities (18). However, these biomarkers have limited sensitivity; thus, by themselves, have not proven to be accurate enough to replace cystoscopy or even VUC. The idea that the presence or absence of a single molecular marker will aid diagnostic or prognostic evaluation has not proved to be the case. This makes sense when one realizes the complex interactions between various molecules within a single pathway, the cross-talk between molecular pathways, the redundancy of some pathways, and the oligoclonality of many tumors. Thus, an evolution to a more global assessment of cancer and the derivation of multiplex molecular classifier signatures holds more promise.
This study is part of a long-term program to develop robust, multiplex molecular assays for noninvasive diagnosis of bladder cancer. We have previously reported in a pilot study the feasibility of profiling the entire transcriptome of shed urothelial cells (3). In this study, we expanded the genome-wide transcriptome profiling to 92 urothelial samples and identified candidate biomarkers associated with the presence of bladder cancer. A panel of differentially expressed genes was selected from the profile data for further investigation in an independent cohort of 81 subjects. The overall strategy is to profile as many genes in as many samples as is technically and economically feasible, followed by more targeted monitoring of a smaller subset of genes in additional samples using a more accurate and cost-effective technique. This approach has enabled us to identify candidate biomarkers and to refine potential multiplex signatures for future validation studies. To date, we have kept both the controls and bladder cancer cohorts relatively separate and homogeneous. This enables identification of promising biomarkers that can be subsequently tested in more complex sets of samples. Conversely, although more typical with respect to clinical presentation, profiling a very heterogeneous cohort (history of bladder cancer, UTI, gross hematuria, etc.) with a high-dimensional array and that does not have adequate representation of each subtype will compromise the ability to identify statistically significant associations. The next phase in our program will be to investigate the significance of the top-ranked panel of biomarkers identified in this study in larger and more diverse cohorts that have a broader representation of bladder cancer subtypes and controls. From that analysis, an optimal signature and a computational method for deriving a diagnostic score will be locked-down for testing in multisite trials.
The rationale for analyzing the shed urothelial component of urine in this study is 2-fold. The first reason is the fact that we aim to develop noninvasive assays for bladder cancer detection, and reasoned that the analysis of the very component that will be the analyte of a future assay is optimal. However, it is equally important that this analyte enables comparison of samples collected from subjects with nonmalignant conditions. Solid tissue samples are available from surgically excised material, but truly normal tissues are almost never available because it would be unethical to biopsy healthy tissue. The analysis of bladder tissue is also restricted by the occurrence of the "field effect", which describes the phenomenon of widespread urothelial change when a malignant lesion is present. Thus, even apparently normal tissue adjacent to excised tumor is not useful for comparison. This is why tissue-based studies have yielded more prognostic than diagnostic biomarkers of promise to date (19, 20). The analysis of urothelia overcomes this limitation because all individuals shed an appreciable number of cells routinely.
The gold standard for the noninvasive diagnosis of bladder cancer is VUC. VUC relies on the microscopic visualization of shed cancer cells in voided urine. Low-grade tumors and low-stage tumors (Ta–T1) shed less cancer cells into the urine, and therefore, the sensitivity for detecting these tumors by VUC typically ranges from 20% to 40% (21). In the current study, VUC detected 31% and 27% of the subjects with bladder cancer in cohort 1 and cohort 2, respectively. The detection rate of bladder cancer by VUC in our study is in line with its detection rate in our previous studies (3, 22). The amplification of shed urothelial cell molecular content, in this case specific transcripts, by PCR obviously provides a huge advantage over a visual analysis of scant cellular material, but it also enables the detection of informative molecular changes that may not manifest as changes in nuclear morphology—the criteria used for VUC evaluation. Of the genes in our best performing signatures (Table 3), some have been implicated in bladder cancer previously (MMP12, CA9, ERBB2, MMP9, and VEGFA), whereas have been associated with other cancers but not bladder cancer (CCL18, DMBT1, RAB1A, SEMA3D, DSC2, AGT, and CRH). To our knowledge, TMEM45A, SYNGR1, and MXRA8 have not been linked to any cancers previously. Of those implicated in bladder cancer, only ERBB2, CA9, and MMP9 have been investigated at all as candidate biomarkers in urine-based assay for cancer detection, and in this context, only 1 or 2 articles for each marker are available in the literature.
The ability to carry out global gene expression profiling on the minimal material present in shed urothelial cell (3) greatly facilitates the identification and development of potential biomarkers for the detection of bladder cancer in noninvasively obtained patient samples. Urine-based assays have several advantages over tissue as a clinical sample. The noninvasive sampling not only benefits the patient by avoiding discomfort, but also enables repeat and copious sampling. Biomarkers can be monitored at the DNA, RNA, or protein level in the cellular material or the supernatant (18, 23, 24), and although a stand-alone test is desirable, information can also be used to increase the accurate detection power of cytology, cystoscopy, or biopsy evaluation. A European team has also used the urothelial sample approach for molecular classification and signature derivation using a similar strategy (25). In parallel with our microarray profiling, the initial analysis by Mengual and colleagues used bladder washings (barbotage). However, we carried out genome-wide profiling (>47,000 transcripts using Affymetrix U133 Plus 2.0 array) on 92 barbotage specimens to identify promising targets for investigation with a target-specific technique, Mengual and colleagues assessed less targets (384 genes) using quantitative RT-PCR. The 384 genes were selected from a previous study that profiled only 12 solid bladder tissue samples (26) using a cDNA microarray platform. Their analysis identified a 12-gene signature that achieved high accuracy in diagnosing bladder cancer as shown in an independent cohort using quantitative RT-PCR. The 12-gene signature achieved diagnostic accuracy of 89% sensitivity and 95% specificity. Despite significant differences between the current study and the study by Mengual and colleagues, both groups were able to derive molecular signatures that can accurately identify bladder cancer. This shows that a multiplex quantitative PCR test on voided urine sample holds promise as a noninvasive urine-based assay in the evaluation of patients being investigated for bladder cancer. Given the different approaches used to get to a biomarker panel that works in actual voided urine samples, it is not surprising that there is not so much overlap in the diagnostic signature between the 2 urothelia studies. Overlap between gene lists of microarray studies conducted in bladder tissues have also been reported to be less than might be expected (27); however, a number of genes present in our top-ranked lists and models are components of signatures under investigation by others (25). It may be that there is considerable overlap of signature genes at a higher level, either as components of common signaling pathways or through a commonality of function.
We recognize that our study has several limitations. First, as a tertiary-care facility, we tend to see more high-grade, high-stage disease, which is reflected in our study cohort. To further confirm the robustness of our multiplex signature, subsequent studies must assess larger cohorts that include subjects with low-grade, low-stage disease. We also did not have complete smoking data for all subjects in the cohort and, therefore, an association with smoking history was not possible. Third, processed, banked urines were analyzed. Urines were centrifuged and separated into cellular pellet and supernatant before storage at −80°C. It is feasible that freshly voided urine samples may provide different results. We are currently investigating the performance of selected biomarkers in urines processed via a number of different protocols, including freshly voided urines. The use of Affymetrix GeneChip arrays and TLDA avoids a lot of typical technical errors, but in future studies we will also include more technical replicates to assess intrasample variation. Next, the sensitivity of VUC in our cohort of predominantly high-grade (grade 3) disease (31% and 27%) was lower than would be expected. This calls into question the known interobserver variability of interpreting VUC. In subsequent studies, we will use 2 cytopathologists to interpret these results. Finally, our second cohort of 81 cases is relatively small and the 2 groups that comprised the 81 subjects were relatively homogeneous, that is, either active cancer or control cases with no active cancer, no history of cancer, no urinary tract infection, no urolithiasis, and no gross hematuria. Thus, we were not able to assess sensitivity/specificity of our biomarkers among different stages/grades.
In this discovery phase study, we have identified biomarkers and diagnostic signatures that achieved higher values of sensitivity and specificity than previously reported for any current or proposed test, confirming the feasibility of developing diagnostic assays that use shed urothelia as analyte. The development of an accurate, noninvasive bladder cancer detection assay would clearly benefit both the patient and the health system. If the assay is accurate and reliable enough, then it would be applicable for not only diagnosing and surveilling the bladder cancer, but also for screening at-risk asymptomatic individuals. Larger, confirmatory studies are now under way to further refine and validate biomarker panels before the development of clinically applicable urine-based assays.
Disclosure of Potential Conflicts of Interest
V. Urquidi is employed in NonaGen Bioscience Corporation as a CSO. S. Goodison has ownership interest (including patents) in provisional patent as a co-inventor. C.J. Rosser has ownership interest (including patents) in Nonagen Bioscience Corporation. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: V. Urquidi, S. Goodison, C.J. Rosser
Development of methodology: S. Goodison, C.J. Rosser
Acquisition of data: V. Urquidi
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): V. Urquidi, S. Goodison, Y. Cai, Y. Sun, C.J. Rosser
Writing, review, and/or revision of the manuscript: V. Urquidi, S. Goodison, Y. Cai, Y. Sun, C.J. Rosser
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S. Goodison, C.J. Rosser
Study supervision: V. Urquidi, S. Goodison, C.J. Rosser
Acknowledgments
The authors thank the 173 subjects that comprised the study cohorts of this study.
Grant Support
Financial support for this work was provided, in part, by the National Cancer Institute under grant R01CA116161 (to S. Goodison), the Flight Attendant Medical Research Institute grant (to C.J. Rosser), and the Florida Department of Health JEK Grant 1KT-01 (to C.J. Rosser and S. Goodison).
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.