Gene Expression Profiling in Cervical Cancer: An Exploration of Intratumor Heterogeneity

Materials and Methods

Specimen collection. Sixty-six flash-frozen punch biopsies were obtained from 16 patients with cervical cancer who were undergoing examination under general anesthesia as part of their pretreatment evaluation for cervical cancer. From each patient, four to seven biopsies were obtained from different regions of the tumor. Immediately upon removal from the patient, the specimens were bisected; one-half was placed in a storage medium (optimal cutting temperature compound) for conventional histopathologic examination, then immediately flash-frozen in liquid nitrogen. The other half of each sample was flash-frozen without any storage medium. The specimens were subsequently evaluated by a pathologist (J. Schwock) based on the morphology of the frozen H&E-stained tissue sections cut from the optimal cutting temperature–embedded material. The total cell content—including stroma and tumor cells—was estimated for all 66 tissue samples using a light microscope. Only samples which contained >50% tumor cells were considered for further analysis. According to this, 33 tissue samples from 11 different patients were included with a range of tumor cell fractions of 60% to 100% (Supplementary Table S1; Supplementary Fig. S1).

Pretreatment evaluation for these patients included bimanual rectovaginal palpation of the tumor, intratumor oxygenation, and interstitial fluid pressure measurements as previously described (24, 25). This study was approved by the Research Ethics Board of the University Health Network, and written consent was obtained from all participating patients.

RNA purification. Total RNA was purified using the Absolutely RNA Microprep Kit (Stratagene, La Jolla, CA). Briefly, tumor tissue was lysed in guanidine thiocyanate, mixed with ethanol, and RNA was captured on a silica-based fiber matrix within a microspin cup. RNA was then washed with saline buffers, DNA contamination was removed with DNase, and RNA was then eluted in RNase-free water. The quality of the purified RNA was assessed by analyzing 200 pg of each sample using a Bioanalyzer 2100 (Agilent, Palo Alto, CA). Only samples with a 28S/18S ribosomal peak ratio of 1.8 to 2.0 were considered suitable for labeling.

Microarray hybridization. The Human Genome U133A Plus 2.0 Gene Chip (Affymetrix, Santa Clara, CA), which contains 54,675 ProbeSets, representing 24,325 distinct UniGene clusters, was used for this study. A total of 1.5 μg of purified total RNA template was reverse-transcribed to generate double-stranded cDNA. Following second-strand cDNA synthesis, biotin-labeled antisense cRNA was generated by in vitro transcription. Next, 15 μg of each generated cRNA preparation was fragmented and hybridized to an oligonucleotide array. Automated washing, staining, and scanning was done according to the manufacturer's protocols.

Normalization and analysis. The raw data were preprocessed using the GCRMA algorithm (26). First, the expression signals of the perfect match probes were corrected for optical noise and nonspecific binding by incorporating mismatch probe information. Next, individual probe intensities were smoothed through quantile normalization (27). Finally, expression values for each ProbeSet were generated using a median polish (28). This algorithm was implemented in the GCRMA package (v1.1.3) of the Bioconductor open source library (29) for R (version 2.0.1). Raw and normalized microarray data have been deposited in the GEO repository at NCBI under accession GSE5787.

Filtering process and clustering. The resulting expression values were filtered to remove low-intensity signals from unexpressed genes, which represents experimental noise. Because only samples from female tumors were analyzed, signals of Y chromosome genes should reflect nonspecific hybridization. Therefore, Y chromosome genes were used to estimate a threshold that defines whether or not a gene is expressed. Accordingly, 34,675 ProbeSets with GCRMA-normalized signal intensities <4.0 were excluded (Supplementary Fig. S2). Unless noted otherwise, all additional analyses employed only the 20,000 ProbeSets with normalized intensities >4.0. Unsupervised hierarchical agglomerative clustering was done using the cluster package in R software (version 2.0.1) and the Euclidian distance as a measure of dissimilarity.

ANOVA. A variance-component analysis of the expression values for the 20,000 putatively expressed ProbeSets was done. The variance within a patient (W = variance due to differences within a tumor), the variance between patients (B = variance due to differences between patients) and the total variance (T = W + B) were produced. The ratio W/T was then calculated and used as a measure of intratumor heterogeneity.

Associations between measures of intratumor and intertissue variability. To determine if the intratumor heterogeneity of a gene's expression might be related to its regulatory heterogeneity, we attempted to correlate the W/T values determined in our study with estimates of intertissue variability obtained from the GNF Gene Atlas (30). Three measures of variability were considered: W = within-sample variance; B = between-sample variability; T = total variance (within-sample + between-sample variance); W/T, a measure of intrasample heterogeneity, and B/T, a measure of intertissue variability. In addition, we also considered the mean intertissue expression level from the GNF Atlas and a Spearman correlation analysis was done between these measures.

Functional gene ontology analysis. Ontological analysis, using the GOMiner software (build 148; ref. 31), was employed to understand the functional relevance of intratumor heterogeneity. First, the set of 20,000 ProbeSets was ranked according to their W/T ratio. Next, this ranked list was broken into 10 deciles, each containing 2,000 ProbeSets. Employing the NetAffx annotation for each ProbeSet (annotation date: July 25, 2005), the distinct genes present in each decile were identified. GOMiner was used to determine the probability of functional enrichment in each decile relative to the entire set of 20,000 ProbeSets. The identified categories were grouped manually and color-coded according to the percentage of all annotated ProbeSets in any given decile.

Real-time quantitative PCR analysis. Real-time PCR (RT-PCR) amplification was done on 30 of the 33 samples, using primer sets for ACTB (as a normalization gene), DUSP1, CD55, IRAK1, CSTA, IL8, HIF-1α, and VEGF genes. Three samples were excluded due to insufficient material. β-Actin was used as a control for the amounts of cDNA generated from each sample. Synthesis of the first-strand cDNA was carried out using SuperScript First-Strand Synthesis System for QRT-PCR (Invitrogen, Carlsbad, CA). The RT product (1-3 μL) was then amplified for 40 cycles (2 minutes at 50°C, 15 minutes at 95°C, 15 seconds at 94°C, 30 seconds at 58°C) followed by an extension of 30 seconds at 72°C. Each assay was repeated thrice, and the mean ΔC_T values were used for further calculations. The sequences of the PCR primer pairs were aligned to the mRNA, and to the Affymetrix ProbeSet for these genes (Supplementary Table S2). A Spearman correlation analysis for signal intensity of the microarray data and ΔC_T values of quantitative RT (QRT)-PCR data was then done.

Estimation of the number of biopsies required as a function of gene expression heterogeneity. The total variance was divided between variance due to patient heterogeneity (σ̂_B²) and variance due to tumor heterogeneity (σ̂_W²). When more than one sample per patient is analyzed, the variance of the mean value per patient decreases as the number of replicates per patient increases. Thus, the variance of the mean per tumor when k samples are analyzed is: . The total variance in this case is: .

Results

Details of patient and tumor characteristics are provided in Table 1 , describing features typical of cervical cancer patients who are treated with chemoradiotherapy. On average, three biopsies were obtained from different regions of the tumor from each patient.

Variance analysis. The ratio W/T, which is a measure of the intratumor heterogeneity in gene expression, was calculated for each of the 20,000 putatively expressed ProbeSets. Approximately equal numbers of ProbeSets displayed low intratumor heterogeneity [W/T < 0.1; 448 ProbeSets (2.2%)] and high intratumor heterogeneity [W/T > 0.9; 537 ProbeSets (2.7%)]. The majority (95%) of ProbeSets fell between these two extremes, but with a global tendency towards lower intratumor heterogeneity. In total 14,473 ProbeSets (72.4%) had higher intertumor than intratumor heterogeneities (W/T < 0.5; Fig. 1 ); thus, for most genes, there was a greater variation between patients' tumors than within any single individual's tumor.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

W/T from 20,000 ProbeSets. The distribution of the 20,000 ProbeSets as a function of W/T values for all 11 patients' tumors. An increasing value of W/T represents an increasing degree of intratumoral heterogeneity in gene expression within patients' cervical cancer specimens.

In order to explore whether the intratumor heterogeneity is dependent on the expression level of a gene, we examined the relationship between ProbeSet signal intensity and W/T. The W/T values were dichotomized according to the mean ProbeSet signals. A significant trend was observed with higher mean signal intensities being represented by genes with lower W/T ratios (P < 0.0001, χ² test; Supplementary Table S3), suggesting that genes with higher expression levels are associated with reduced intratumor heterogeneity. However, it must be stressed that the very small P value is mainly due to the very large number of genes analyzed.

Validation with QRT-PCR. In order to internally validate the results obtained from the microarray data, we selected seven representative genes with different W/T ratios for further evaluation (Table 2 ). Three genes were selected which had low intratumor heterogeneity (low W/T values), included DAF, CSTA, and IRAK1. One gene with high intratumor heterogeneity (high W/T ratio) was chosen (DUSP1), and three genes functionally linked to hypoxia and angiogenesis were also selected (HIF-1α, IL8, and VEGF). The variation of each of these genes within and between patients is illustrated graphically in Fig. 2 , indicating that genes with low intratumor heterogeneity (low W/T ratios) indeed do cluster closely within each patient's tumor, whereas DUSP1 (with high W/T) has different degrees of expression intensities within many patients' tumors.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Scatter plots demonstrating intratumoral and intertumoral variability for selected genes. Scatter plots of signal intensities of ProbeSets representing genes with W/T < 0.15 (DAF, CSTA, and IRAK1), with W/T >0.75 (DUSP1) and genes related to angiogenesis (HIF-1α, IL8, and VEGF). The x-axis denotes the individual patients (n = 11); each dot indicates the expression intensity of that gene, from each sample obtained from the corresponding patient (n = 33). The y-axis represents normalized signal intensity. Dots that cluster together show minimal intratumor variation (e.g., DAF); dots that are further apart show a higher degree of intratumor variation (e.g., DUSP1).

Table 2 provides further details on each of these genes, including their putative function, and chromosomal locations. The W/T ratios are also provided for both the microarray signal intensity, and the QRT-PCR ΔC_T values, along with correlations between the microarray and QRT-PCR measurements in terms of the W/T values. A strong correlation (correlation coefficient > 0.7) was observed for DAF, IL8, and VEGF; satisfactory correlations (0.5-0.7) were observed for CSTA, IRAK1, and DUSP1. No correlation was observed, however, for HIF-1α.

Two samples were identified as “outliers” in all QRT-PCR results (Fig. 3 ). When these samples were re-examined in greater detail, we went back to the original QRT-PCR data, samples labeled “12.5” and “4.1” gave unreliable results, indicating that RNA may have been damaged after repeated thawing and freezing. After excluding those two samples, the correlation coefficients were 0.828 for DAF, 0.834 for IRAK1, 0.770 for VEGF, 0.717 for DUSP1, 0.822 for IL8, and 0.669 for CSTA, indicating a very good correlation between microarray and QRT-PCR results (Table 2; Fig. 3).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Correlation between microarray and QRT-PCR results. Histograms of the correlation between signal intensity from the microarray analyses, with the QRT-PCR results for DAF, CSTA, IRAK, HIF-1α, IL8, VEGF, and DUSP1. The x-axis represents the ΔC_T values obtained from QRT-PCR, the y-axis represents the normalized ProbeSet intensities obtained from the microarray experiments.

Hierarchical clustering. An unsupervised clustering was done to determine the extent of similarity in global gene expression pattern for each patient's tumor, and whether different patients are distinguishable on the basis of this profile. When a cluster analysis was done on all 33 samples, all but one specimen (s10_7) clustered perfectly within each corresponding patient, indicating that the expression pattern is very similar within each patient (Fig. 4A ). This clustering suggests that there is a unique expression profile that characterizes each patient's tumor. Furthermore, when the cluster analysis was repeated, but restricted to ProbeSets with low intratumoral heterogeneity (W/T < 0.5), this resulted in perfect clustering of samples within each corresponding patient (Fig. 4B).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Dendrogram of 33 samples clustered by gene expression patterns. The 20,000 ProbeSets from the 33 samples obtained from 11 cervical cancer patients, with mean signal intensities ≥4 (A); the same analysis, but including only ProbeSets with W/T values of <0.5 (B). Each branch of the tree represents an individual sample from each patient (samples are labeled: “s”, patient ID, underscore, and sample number for that patient). All samples, except one, clustered perfectly with the corresponding patient, indicating greater similarity in gene expression pattern within patients, than between different patients (A). This is even more clearly shown for ProbeSets whose expression heterogeneity is limited to W/T values of <0.5 (B).

Prediction of expression heterogeneity. In order to determine if the intratumoral heterogeneity in gene expression might be related to its regulatory heterogeneity, we attempted to correlate the W/T values determined in our study with estimates of intertissue variability obtained from the GNF Gene Atlas (30). We selected only those ProbeSets present identically on the array platforms used for both studies (12,533 of the 20,000 putatively expressed ProbeSets). The Spearman correlation coefficients between these measures are shown in Supplementary Table S4. Only very weak correlations were observed, suggesting that intertissue heterogeneity does not account for the variations in gene expression intensity within each tumor.

Functional analysis of low versus high W/T genes. The 20,000 ProbeSets in our data set represent 11,141 distinct genes, based on the Affymetrix annotation. We wished to determine if intratumor gene expression heterogeneity might be related to distinct biological functions, hence, we tested the relationship between W/T values as a function of Gene Ontology (GO) annotations, using the GOMiner software. A total of 61 GO categories seem to be enriched for genes with very high intratumor homogeneity (P < 0.001 in the either of the two lowest deciles), or with high intertumor heterogeneity (P < 0.001 in either of the two highest deciles). These categories were grouped manually according to functionally related groups (see Fig. 5 ).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Ontological analysis of genes with high W/T (top 10%) and genes with low W/T (bottom 10%). To determine if intratumor variability was related to gene function, we employed the GOMiner program for ontological analysis. A subset of categories which showed enrichment in highly homogeneous or heterogeneous genes at the P < 0.001 level. Color scale coding reflects the percentage of genes with a given annotation found in a specific decile. For example, fully 25% of ProbeSets with the GO annotation “MHC class II receptor activity” had W/T ratios <0.1 of the dataset, whereas none had W/T ratios >0.9. By contrast, >15% of ProbeSets annotated as “RNA splicing” had W/T >0.9, whereas only 4.4% had W/T ratios >0.9.

Transcripts with low W/T (<10%) were characterized by a statistically significant enrichment of genes with immune-related function, i.e., antigen processing and MHC class II receptor function (P < 0.0001; Fig. 5). Genes with high W/T were characterized by a significant tendency towards nuclear localization and nucleic acid binding (both P < 0.0001). Fully 29% (85 of 293) of genes with W/T > 0.95 were involved in nuclear localization, in comparison with only 19% (2,148 of 11,141) of genes in the complete data set (Fig. 5).

Estimation of the number of biopsies required as a function of gene expression heterogeneity. A model to estimate the necessary number of biopsies required to obtain reliable expression results for each gene based on the W/T of that gene, was developed (Fig. 6 ; refs. 32, 33). Using this model, the number of biopsies and the corresponding reduction in W/T were estimated. A stringent measure to reduce heterogeneity would be W/T ≤ 0.15, shown as black circles in Fig. 6 (32, 33). The model predicts that as the W/T value for a gene increases (greater intratumoral heterogeneity), a larger number of biopsies would be necessary to obtain a representative “whole tumor” expression for that specific gene. Hence, for genes with a W/T ≤ 0.15, a single biopsy would suffice; for genes with W/T ranging from 0.16 to 0.25, two biopsies are necessary. For genes with W/T values between 0.26 and 0.30, three biopsies would be required; for genes with W/T 0.3 to 0.4, four biopsies would be necessary. Finally, for highly heterogeneously expressed genes with W/T values >0.5, six or more biopsies would be required.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

A model that provides an estimate for the number of biopsies necessary to obtain a reliable read-out for each gene, as a function of its W/T value. The y-axis represents the number of biopsies required; the x-axis represents the W/T for each gene. •, number of biopsies that would be necessary in order to reduce heterogeneity to ≤0.15, in order to obtain a more reliable estimate for gene expression. For genes with a W/T ≤ 0.15, a single biopsy would suffice; for highly heterogeneous genes with W/T values >0.5, six or more biopsies would be required.

Discussion

The advent of high-throughput technologies has enabled the definition of patients' tumors as a function of their gene expression profile (34). One concerning aspect of the literature is that different studies report different sets of potentially prognosticating genes, with little overlap among these reports (7, 8). There are many reasons for this variation, but one possibility that has been understudied is the issue of gene expression heterogeneity within patients' tumors.

Previous studies have shown significant intratumor heterogeneity in human cervical cancer, with variations in chromosomal aberrations, microvessel density, pO₂, interstitial fluid pressure and protein expression in different regions within a single patient's tumor (11–19). Based on these many reports, our hypothesis was that there would be significant intratumoral heterogeneity in gene expression profiles in patients with locally advanced cervical cancers. Our data, in fact, shows that the majority of genes were expressed relatively homogeneously, in that 72% of the ProbeSets had intratumoral variations of <0.5 of the total variance (Fig. 1). In an unsupervised cluster analysis, all samples except one grouped perfectly with their corresponding patient (Fig. 4), again indicating that most of the variability in fact occurred in-between different patients' tumors.

There remained, however, significant heterogeneity among a subset of genes. An ontological analysis revealed that genes whose expression varied significantly within a tumor (high W/T), tended to fall into two broad categories; those involved with transcriptional regulation (at the level of expression as well as RNA splicing), and those involved with cellular metabolism. The high intratumor variability of these genes may reflect the changing activity of different cells within the tumor. In this context, it is worth noting that microarray data presents only a snapshot of cell activity, which obviously could change over time.

In contrast to the intratumor variability in the expression of transcription-related genes, translation-related genes, and in particular, the ribosomal genes displayed consistent expression throughout a tumor and were among the least variable genes. In addition, genes involved in cell-mediated immune response, including antigen processing and presentation, and MHC class II receptor activity genes, were also homogeneously expressed throughout the tumor. This is intriguing, as cervical cancer has a specific viral etiology, and cellular immunity likely plays a major role in cervical cancer carcinogenesis (35). In our study, genes with MHC class II receptor activity, as well as genes involved in antigen processing and presentation, were found to be one of the most homogeneously expressed genes, which could be a reflection of host factors, and therefore, less susceptible to regional differences within a tumor.

In an effort to derive some understanding of the implications of our data, we provide an estimate for the number of biopsies that would be necessary, in order to minimize the effect of intratumoral gene expression variability within a single patient's tumor. It has been previously suggested that a reliability coefficient of at least 0.85 would be adequate to render clinical decisions based on results of diagnostic tests (33). Applying this standard to our data, this would imply that only genes with a W/T < 0.15 (n = 1,536) could be reliably estimated from a single biopsy. Consequently, genes with higher W/T (n = 18,464) would require a larger number of biopsies in order to provide a reasonable estimate of their “true” level of gene expression (see Fig. 6). For example, hypoxia-inducible factor 1α (HIF-1α) had a W/T of 0.19; hence, two biopsies from the same patients' tumor would be necessary to obtain a more reliable estimate for this gene. In contrast, VEGF with W/T = 0.70 would require >10 biopsies. If a less stringent estimate of W/T is required (recognizing that for many useful prognostic factors, W/T may range up to 0.5), the model presented in Fig. 6 could be used to estimate the number of biopsies required, based on the initial and required W/T. However, an optimum W/T cannot be chosen on statistical reasoning alone. In the clinical world, the feasibility of taking many biopsies from one patient is a restrictive factor.

The variable proportion of tumor cells in each sample may contribute to heterogeneity in mRNA expression levels. We therefore chose to assay microarray samples with a tumor fraction of at least 50% (Supplementary Table S1; Supplementary Fig. S1). The distribution of tumor fractions is very narrow, with 31 of 33 samples having tumor fractions in the interval (0.6, 0.9). This narrowness should greatly help to reduce the influence of this confounding variable on our analysis of tumor variability, but also limits our ability to distinguish the effects of tumor heterogeneity from those of tumor cell fractions. In our data, no clear relationship between expression signal and tumor cell fraction could be observed (Supplementary Fig. S3).

Tumor sample size may also contribute to intratumor heterogeneity and a small needle biopsy may lead to a less representative measure of tumor gene expression compared with a much larger sample size. In our study, all samples were punch biopsies from similar sizes, although the exact volumes of the samples were not available.

Technical variability may have a negligible effect on the intratumor heterogeneity in our study because Affymetrix arrays are extremely highly correlated (>99%). Extensive data on this fact was presented with the “Latin Square” spike-in experiment by Affymetrix (36).

Clearly, our data have interesting implications for the continuing development of gene expression profiling as providing prognostic or predictive information for human cancers (37). A useful prognostic marker would not only have to correlate with patient outcome, but it also ideally should be expressed homogeneously within any single patient's cancer, so that sampling “error” would not become a confounding variable. This implies that genes with low intratumor heterogeneity (low W/T values) would be preferentially more useful as biomarkers than those with higher intratumor heterogeneity (higher W/T values). The incorporation of this additional W/T variable into the algorithms for the selection of prognostic subsets would be an important yet challenging computational issue.

The observations from this study can only be considered as preliminary because we studied a relatively small number of patients; our results need to be confirmed by an independent cohort of patients with cervical cancer. The issue of intratumoral heterogeneity has been evaluated recently for gastric cancer (38), whereby very little intratumoral heterogeneity in gene expression was observed. However, this study was even more limited in that only six patients were evaluated, and only two samples from each patient's tumor were subjected to the Affymetrix U133A oligonucleotide analysis. In contrast, another study reported significant heterogeneity in the expression of stress and hypoxia-activated genes in lung cancer surgical specimens (39), which corroborates our own observation that genes involved in transcription or metabolism are expressed variably in different regions of the tumor.

The finding of increased heterogeneity among the low abundant transcripts is in agreement with previous work by O'Sullivan et al. (40), who also studied intratumor heterogeneity in gene expression by microarray techniques. They claim that technical factors could not account for all the increased variability of these genes, and that biological factors were also probably important.

In summary, the majority of genes are expressed relatively uniformly within any single patient's cervical tumor. However, a subset of genes involved in metabolism and transcriptional regulation can be expressed quite variably within a single cancer, indicating that for such genes, multiple samplings would be required in order to account for this heterogeneity, to fully understand their true prognostic or predictive value for cervical cancer.