An ethnically diverse cohort of patients was studied using samples collected from the University of Utah Health Sciences Center, from the University of North Carolina, from Thomas Jefferson University, from the Maine Medical Center, and from the University of Chicago. Patients provided written acknowledgement of informed consent in accordance with institutional and federal guidelines. Samples collected prospectively for microarray and qRT-PCR analyses included 117 invasive breast cancers, one fibroadenoma, five 'normal' samples (from reduction mammoplasty), and three cell lines. Patients were treated in accordance with the standard of care dictated by their disease stage, ER status, and HER2 status. Patient outcome information was collected for up to 118 months (median 21.5 months). The clinical data for the qRT-PCR samples are presented in Additional file 2 (Supplemental Table 1). Publicly available data sets containing 337 samples with long-term follow-up (median 86.7 months) were used to further validate the prognostic significance of the proliferation meta-gene within the context of intrinsic subtypes [7
Sample preparation and first-strand synthesis for qRT-PCR
Nucleic acids were extracted from fresh frozen tissue using the RNeasy Midi Kit (Qiagen Inc., Valencia, CA, USA). The quality of RNA was assessed using the Agilent 2100 Bioanalyzer with the RNA 6000 Nano LabChip Kit (Agilent Technologies, Palo Alto, CA, USA). All samples used had discernable 18S and 28S ribosomal peaks. First-strand cDNA was synthesized from approximately 1.5 μg total RNA using 500 ng Oligo(dT)12–18 and Superscript III reverse transcriptase (1st Strand Kit; Invitrogen, Carlsbad, CA, USA). The reaction was held at 42°C for 50 minutes followed by a 15-minute step at 70°C. The cDNA was washed on a QIAquick PCR purification column and was stored at -80°C in 25 mM Tris, 1 mM ethylenediamine tetraacetic acid at a concentration of 5 ng/μl (concentration estimated from the starting RNA concentration used in the reverse transcription).
Genbank sequences were downloaded from Evidence viewer (NCBI website) into the Lightcycler Probe Design Software (Roche Applied Science, Indianapolis, IN, USA). All primer sets were designed to have a Tm of approximately 60°C, to have a GC content of approximately 50%, and to generate a PCR amplicon <200 bps. Finally, BLAT and BLAST searches were performed on primer pair sequences using the UCSC Genome Bioinformatics database and the NCBI database to check for uniqueness. Primer sets and identifiers are provided in Additional file 2 (Supplemental Table 2).
For PCR each 20 μl reaction included 1 × PCR buffer with 3 mM MgCl2 (Idaho Technology Inc., Salt Lake City, UT, USA), 0.2 mM each of dATP, dCTP, and dGTP, 0.1 mM dTTP, 0.3 mM dUTP (Roche Applied Science), 10 ng cDNA and 1 U Platinum Taq (Invitrogen). The dsDNA dye SYBR Green I (Molecular Probes, Eugene, OR, USA) was used for all quantification (1/50,000 final). PCR amplifications were performed on the Lightcycler software using an initial denaturation step (94°C, 90 seconds) followed by 50 cycles: denaturation (94°C, 3 seconds), annealing (58°C, 5 seconds with 20°C/s transition), and extension (72°C, 6 seconds with 2°C/sec transition). Fluorescence (530 nm) from the dsDNA dye SYBR Green I was acquired for each cycle after the extension step. The specificity of the PCR was determined by postamplification melting curve analysis. Reactions were automatically cooled to 60°C at a rate of 3°C/s and slowly heated at 0.1°C/s to 95°C while continuously monitoring the fluorescence.
Relative quantification by real-time qRT-PCR
Quantification was performed using the LightCycler 4.0 software. The crossing threshold for each reaction was determined using the second-derivative maximum method [15
]. The relative copy number was calculated using an external calibration curve to correct for PCR efficiency and a within-run calibrator to correct for the variability between runs. The calibrator is made from four equal parts of RNA from three cell lines (MCF7, SKBR3, and ME16C) and Universal Human Reference RNA (catalogue number #740000; Stratagene, La Jolla, CA, USA).
Differences in cDNA input were corrected by dividing the target copy number by the arithmetic mean of the copy number for three housekeeper genes (MRPL19, PSMC4, and PUM1) [17
]. After adjusting copy numbers to the reference sample (calibrator) in LCS4, relative copy numbers were imported into a relational database where the data were normalized to the housekeeper genes and were log2
-transformed for further analyses.
Hierarchical clustering was carried out in Cluster analysis using Spearman correlation, median centering by gene and array, and average linkage association [18
]. The clustering was visualized using Treeview. The real-time qRT-PCR relative copy number data for all genes (53 classifier genes and three housekeeper genes) can be found in Additional file 2 (Supplemental Table 3).
Histological assessment of grade was performed for the invasive ductal adenocarcinomas using the Scarff-Bloom-Richardson system. Nuclear grading was determined for tumors in which tubular differentiation could not be assessed (for example, invasive lobular carcinomas). Samples were scored for protein expression at the time of diagnosis and using standard operating procedures established at each institution. Greater than 20% positive staining nuclei was considered positive for the ER and the progesterone receptor. Staining and scoring criteria for HER2 were carried out according to the HercepTest (Dako, Carpinteria, CA, USA).
The 126 samples used for qRT-PCR were also analyzed by DNA microarray (Agilent Human A1, Agilent Human A2, and custom oligonucleotide). Labeling and hybridization of RNA for microarray analysis were performed using the Agilent low RNA input linear amplification kit, but with one-half of the recommended reagent volumes and using a Qiagen PCR purification kit to clean up the cRNA. Each sample was assayed versus a common reference sample that was a mixture of Human Universal Reference total RNA (Stratagene, La Jolla, CA, USA) enriched with equal amounts of RNA from the MCF7 and ME16C cell lines. Microarray hybridizations were carried out on Agilent Human oligonucleotide microarrays using 2 μg Cy3-labeled 'reference' sample and 2 μg Cy5-labeled 'experimental' sample. Hybridizations were carried out using the Agilent hybridization kit and a Robbins Scientific '22k chamber' hybridization oven (Robbins Scientific, Sunnyvale, CA, USA). The arrays were incubated overnight, washed once in 2 × SSC and 0.0005% Triton X-102 (10 minutes), washed twice in 0.1 × SSC (5 minutes), and were then immersed into Agilent Stabilization and Drying solution for 20 seconds.
All microarrays were scanned using an Axon Scanner 4000A (Axon Instruments, Inc, Foster City, CA, USA). The image files were analyzed with GenePix Pro 4.1 (Axon Instruments) and were uploaded into the UNC Microarray Database at the University of North Carolina at Chapel Hill, where a lowest normalization procedure was performed to adjust the Cy3 and Cy5 channels [19
]. All primary microarray data associated with this study are available at the UNC Microarray Database and have been deposited in the GEO under accession number GSE2607.
Selecting genes for real-time qRT-PCR
We developed a real-time qRT-PCR assay using 53 genes that were selected due to their importance in making 'intrinsic' subtype distinctions and/or their association with cell proliferation (see Additional file 2, Supplemental Table 2). The statistical selection of 'intrinsic' genes involved using 45 before-therapy and after-therapy samples derived from the data set presented in Sorlie and colleagues (see Additional file 2, Supplemental Table 4 for the list of 45 pairs) [5
]. The two-color DNA microarray data were downloaded from the Internet and the R/G ratio (experimental/reference) for each spot was normalized and log2
-transformed. Missing values were imputed using the k-NN imputation algorithm described by Troyanskaya and colleagues [20
Using an 'intrinsic' analysis [3
] we identified 550 microarray elements/spots from the data set presented by Sorlie and colleagues [5
]. We then applied the 'intrinsic' genes to identifying molecular subtypes within a completely independent data set of early-stage breast cancers [7
]. Common elements between the data sets were found after translating the gene annotation from each data set to UniGene Cluster IDs using the SOURCE database [21
]. Following the algorithm outlined by Tibshirani and colleagues [22
], we hierarchical clustered the 97 samples from van 't Veer and colleagues' study [7
] using a common set of 350 genes and assigned intrinsic subtypes (Luminal, HER2+/ER-, Basal-like, or Normal-like) based on the sample-associated dendrogram. Finally, we identified genes that optimally distinguished the four subtype classes using a version of the gene selection method first described by Dudoit and Fridlyand [24
], where the best class distinguishers are identified according to the ratio of between-group to within-group sums of squares. After scoring genes in this manner, 10-fold cross-validation was performed with a nearest centroid classifier, resulting in a list of 41 genes that gave the highest prediction accuracy when compared with the entire set of 350 genes.
We successfully developed qRT-PCR assays for 37 out of the 41 genes identified from the 'intrinsic' analysis (minimal intrinsic list). The genes PGR and EGFR were also included in the qRT-PCR assay, despite not being statistically selected in the 'intrinsic' analysis, because of their value in predicting therapy response and their strong association with ER-positive and ER-negative tumors. Finally, we tested 14 proliferation genes because of their importance in prognosis.
The stability of our hierarchical clustering classifications was tested using a k
-means algorithm implemented in Cluster 3.0 and using Consensus Cluster [25
] implemented in GenePattern. For the k
-means we ran 1,000 trials (K
= 4 and K
= 5) and used the Euclidean distance as the similarity metric. Consensus clustering was also performed for 1,000 runs using a 'subsampling' of 0.8, so that 20% of the samples were left out of each run.
Combining microarray and qRT-PCR datasets
We used distance-weighted discrimination (DWD) to identify and correct systematic biases across the microarray and qRT-PCR datasets [26
]. Prior to DWD, we normalized each dataset by setting the mean to 0 and the variance to 1. After performing DWD, genes in common between the datasets were clustered using Spearman correlation and average linkage association.
Receiver operator curves
In order to determine agreement between protein expression (immunohistochemistry (IHC)) and gene expression (qRT-PCR), a cutoff value for the relative gene copy number was selected by minimizing the sum of the observed false-positive and false-negative errors; that is, minimizing the estimated overall error rate under equal priors for the presence/absence of the protein. The sensitivity and specificity of the resulting classification rule were estimated via bootstrap adjustment for optimism [27
Survival curves were estimated by the Kaplan-Meier method and compared via a log-rank or stratified log-rank test as appropriate. The standard clinical pathological parameters of age (years), node status (positive versus negative), tumor size (cm, a continuous variable), grade (1–3, a continuous covariate), and ER status (positive versus negative) were tested for differences in relapse-free survival (RFS) and overall survival using the Cox proportional hazards regression model. Pairwise log-rank tests were used to test for equality of the hazard functions among the intrinsic classes. Cox regression was used to determine predictors of survival from continuous expression data. All statistical analyses were performed using the R statistical software package (R Foundation for Statistical Computing Vienna, Austria).