Meta-analysis of Colorectal Cancer Gene Expression Profiling Studies Identifies Consistently Reported Candidate Biomarkers

Colorectal cancer, defined as cancerous growths in the colon, rectum, or appendix, is the third most frequent cancer in both males and females in North America (1). This year, an estimated 20,800 Canadians will be diagnosed with colorectal cancer and ∼8,700 will die of it (2). A common area of research interest is the identification of diagnostic biomarkers for early and accurate detection of colorectal cancer (3). Prognostic biomarkers are also being developed to, for example, separate patients who will benefit from adjuvant therapy from those who will not (4) or to determine which patients are at risk for disease recurrence (5). Other studies have focused on understanding cancer progression by identifying differences in gene expression between normal, benign adenoma, and carcinoma stages (6-8).

A common starting point for these studies is the surgical resection of both cancer and normal tissues from patients followed by global expression profiling to determine differentially expressed genes. These studies can result in tens to thousands of such genes, only a small portion of which may actually be of clinical utility. Although an abundance of data comparing the expression profiles of cancerous to normal tissue has been generated, to date, no reliable biomarker has resulted. One explanation for this lack of translational success to the clinic has been the inconsistency in the results of independent studies (1, 9, 10). Explanations for this low overlap include utilization of different tissue resection methods (microdissection, laser capture microdissection, etc.), different expression profiling technologies [cDNA two-channel microarrays, oligonucleotide microarrays, Serial Analysis of Gene Expression (SAGE), etc.], and different analysis methods (multiple correction tests, fold change thresholds, etc.).

Review articles often include lists of genes that have been reported in multiple independent studies. Consistently reported genes are considered to be biologically relevant to colorectal cancer, whereas those reported only sporadically are thought to have resulted from inherent noise or biases in the different platforms and analysis methods employed (11). Although these lists are helpful in summarizing the current biomarker candidates, the statistical significance of the level of overlap is usually not considered. One can imagine randomly choosing genes from each expression profiling platform, randomly labeling them as up-regulated or down-regulated, and observing some overlap due to chance alone. Therefore, it would be useful to determine which differentially expressed genes were consistently reported in independent colorectal cancer expression profiling studies with a statistically significant frequency.

To address these challenges, a meta-analysis method was recently developed by our group and applied to published studies of differentially expressed genes in thyroid cancer (12). Such a meta-analysis ignores the differences between studies, such as the expression profiling platform used, and instead focuses on elucidating consistently reported genes. The meta-analysis method involved a vote-counting strategy in which a gene was ranked according to the number of studies reporting its differential expression, the total number of tissue samples used in the studies, and the average fold change. That study resulted in a panel of 12 differentially expressed genes reported at a frequency highly unlikely to have occurred by chance. The panel contained both well-known thyroid cancer markers as well as some uncharacterized genes, showing the ability of the meta-analysis method to highlight novel candidate biomarkers. With these results in mind, the objective of the current study was to apply the meta-analysis method to colorectal cancer to observe whether a statistically significant level of overlap among studies could be observed and to identify promising biomarkers. Also, we improved the meta-analysis method by dividing genes into semiquantitative categories based on the number of tissue samples to highlight genes that may have shown the greatest fold changes but would have been ranked lower by the original method due to fewer tissue samples studied. We curated published lists of differentially expressed genes from 25 independent studies performing global expression profiling to compare colorectal cancer to normal tissue, adenoma to normal tissue, and colorectal cancer to adenoma tissue. Many genes were consistently reported as differentially expressed in multiple studies and this overlap was highly significant. The list of candidate biomarkers we present here will be a valuable resource to the colorectal cancer research community for further studies.

We queried PubMed for colorectal cancer expression profiling studies published between 2000 and 2007. Only studies using tissue samples obtained from surgical resection of cancerous tumors and/or adenomatous polyps were considered. Studies were divided into three comparison types: cancer versus normal, adenoma versus normal, and cancer versus adenoma. We excluded the limited studies that focused on the microsatellite stability of the tissues, specific Dukes stages, or those comparing cancer to cancer samples to determine prognostic biomarkers. In total, differentially expressed genes from 25 independent studies were collected. Twenty-three studies did expression profiling to compare cancer versus normal tissue samples (Table 1), whereas seven and five studies considered adenoma versus normal (Table 2) and cancer versus adenoma (Appendix 6), respectively.

In the microarray expression profiling studies, differentially expressed genes were represented by an accession ID, HUGO gene name, or Affymetrix probe ID. The sequence identifier was mapped to the National Center for Biotechnology Information Entrez Gene Identifier (Entrez Gene ID; ref. 13) with the aid of custom-developed Perl scripts and the Clone/Gene ID Converter tool (14). For the SAGE study, updated tag to gene mapping data were obtained from SAGE Genie (15).

To estimate the background levels of overlapping studies, we obtained the platform-specific annotation file for each study to identify genes that could potentially be detected as differentially expressed. For commercial platforms, such as Affymetrix and Atlas microarrays, the annotation file was obtained directly from the company Web site. The identifiers in these annotation files were mapped to the corresponding Entrez Gene ID as above to produce a total gene list for each study. Identifiers that could not be mapped to an Entrez Gene ID were ignored. To obtain a total gene list for the SAGE study, all gene names in the tag to gene mapping data from SAGE Genie were mapped to Entrez Gene IDs. For studies that used platforms in which an annotation file could not be obtained, such as the custom cDNA microarrays and some of the oligonucleotide microarrays, an approximation approach was used in which the appropriate number of Entrez Gene IDs was randomly chosen from the combined gene lists from the other platforms.

To determine if the level of overlap among the studies was significant, we did simulations as described previously (12). Briefly, Perl scripts were created to perform Monte Carlo simulations. In each of the 10,000 permutations, the appropriate number of Entrez Gene IDs from the total gene list of each study was randomly chosen and each ID was randomly labeled as “UP” for up-regulated or “DOWN” for down-regulated. We used an “all-or-none” approach in which the level of overlap for a particular gene was only considered if all the independent studies reporting its differential expression agreed on the direction. The level of overlap among studies in each permutation was counted as in the real analysis. On completion of the permutations, a distribution of overlap results from the simulations was determined and a P value was estimated by comparing the overlap from the simulations to the actual level of overlap in the real data. Significance was defined at P < 0.05. Similar to the previous meta-analysis (12), genes were ranked according to three criteria in the following order of importance: (a) level of overlap (that is, listing the same gene as differentially expressed with a consistent direction of change), (b) total number of samples for overlapping studies, and (c) average fold change reported by the studies in agreement. We further subdivided the genes into three categories using a semiquantitative scale: lowest (Q1), moderate (interquartile range), and greatest (values greater than those in Q3) tissue sample sizes. This improvement over the previous published method gives greater importance to the average fold change criteria for ranking genes when total sample numbers are similar.

Of the total 8,176 differentially expressed genes reported in the 25 studies (4,273 up-regulated and 3,903 down-regulated), 7,287 (89.1%) could be mapped to an Entrez Gene ID (3,822 up-regulated and 3,465 down-regulated). In the cancer versus normal and adenoma versus normal comparisons, significant overlap was observed. No such significance was seen in the cancer versus adenoma comparison (Table 3), although each individual study identified differentially expressed genes (refs. 6, 8, 17, 19, 35; see Appendix 6).

We present the results from the cancer versus normal comparison as an example. The simulations showed that the amount of overlap in this comparison was highly significant (P < 0.0001), with 573 genes reported as differentially expressed with consistent direction of change in at least two studies (multistudy genes; Fig. 1). There were 175 multistudy genes that were reported with inconsistent direction of differential expression. Thus, the majority of multistudy genes (76.6%) that were reported as differentially expressed agreed on the direction, even for large numbers of studies.

From the Monte Carlo simulations, an average of 258.30 (95% confidence interval, 258.16-258.45) genes would be expected to have an overlap of 2, whereas the actual data contained 410. An average of 18.37 (95% confidence interval, 18.33-18.42) genes would be expected to have an overlap of 3 compared with 95 in the actual data. For an overlap of 4, the simulation produced 1.14 (95% confidence interval, 1.13-1.15) genes, whereas the actual data contained 30 genes. Overlaps of 5, 6, and 7 were observed in the simulations but with averages of less than one hundredth of a gene. In 10,000 permutations, the simulations never produced an overlap greater than 7, whereas two genes had an overlap of 9 and one gene had an overlap of 11 in the real data. Although the total number of genes with an overlap of 2 was still very significant, we present here only the genes reported by three or more studies, as we deemed these to be the most reliable (Appendices 2-5). Additional information on the results appears in the Appendices (online only).

A logical solution to the problem of lack of agreement among expression profiling studies in colorectal cancer is to determine the overlap among many studies using different platforms and observe which genes are consistently reported as differentially expressed. These genes likely show biological relevance to the tumorigenesis of colorectal cancer, as opposed to sporadically reported genes, which may be false positives.

Meta-analyses have been done previously to determine differentially expressed genes in colorectal cancer (1, 9, 10). However, these studies and others usually do not consider whether the level of overlap observed is statistically significant. In the newest version (3.0) of the cancer profiling database Oncomine (16), a meta-analysis tool was implemented to compare results from independent studies. However, Oncomine presently contains raw data for eight colorectal cancer profiling studies, only two of which would qualify for our study (17, 18), because they were the only studies that performed at least one of the three comparisons of interest. As discussed previously, our meta-analysis method is useful when raw data are unavailable for consistent reanalysis, which is usually the case (12). However, one limitation of our method is that a measure of confidence cannot be assigned at the gene level, such as from calculating a true combined fold change or P value. Thus, in order for more powerful meta-analysis methods to be applied to colorectal cancer profiling studies, researchers should be encouraged to make public their raw data so that they may be included in repositories such as Oncomine.

By applying this method to a near comprehensive collection of colorectal cancer expression profiling studies, we were able to determine the genes that were reported with a statistically significant frequency. As an extension of the meta-analysis method, we categorized some genes according to their total number of tissue samples as lowest (Q1), moderate (interquartile range), and greatest (values greater than those in Q3) instead of using absolute numbers. This allowed the average fold change criterion to have a greater effect on the gene rank in cases where total sample sizes were similar. In the original version of the meta-analysis method, fold change rarely had any effect on rank.

We observed that for the cancer versus normal and adenoma versus normal comparisons, genes were consistently reported as differentially expressed and that the level of overlap was statistically significant. The results from the cancer versus adenoma comparison were not significant, suggesting that the number of multistudy genes in the five studies could have been observed due to chance. Determining the significance of overlap among studies provides another filtering step to remove false-positive genes from further consideration. When ignoring the significance of the observed overlap, one may be misled by multistudy genes. For example, without knowledge of the statistical significance, one may reason that the multistudy genes in the cancer versus adenoma comparison are biologically relevant, although this decision cannot be reasonably made because the observed level of overlap may be due to chance alone.

In the cancer versus normal comparison, a total of 748 genes were reported as differentially expressed in at least two independent studies. Although the majority of these genes were reported as differentially expressed in the same direction, 175 (23.4%) genes were not. Of these 175 genes, 132 (75.4%) were reported in two studies, 32 (18.3%) were reported in three studies, 8 (4.6%) were reported in four studies, 2 (1.1%) were reported in five studies, and 1 (0.6%) was reported in six studies. There are many potential explanations for these observed inconsistencies. Firstly, one limitation with such meta-analyses is the overgeneralization of comparisons. Although every effort was made to ensure that each study included in each of the three comparisons were comparable, there are bound to be inconsistencies due to the lack of relevant clinical data being reported in each of the studies. For example, in the cancer versus normal comparison, SLC26A3 was reported as down-regulated in five studies (8, 17, 19-21) but up-regulated in one study (22). The five studies that reported this gene as down-regulated did not specify the microsatellite status of the colorectal cancer tissue samples being used, whereas the one study that reported the up-regulation of this gene used a mixture of microsatellite stable and unstable tissue samples. Other than microsatellite stability, other clinical features, such as the specific portion of the colon where the tissue samples were taken (9), may affect the direction of differential expression. Thus, due to the lack of these clinical data, it is difficult to determine whether the results of each independent study are truly comparable with each other. Conversely, if these clinical data were more readily available, more specific comparisons, such as microsatellite-stable colorectal tissue samples taken from male patients versus paired normal mucosa, could be done.

A related explanation for why some genes were reported as differentially expressed in an inconsistent direction is the heterogeneity in the tissue samples used. The independent studies experimented on tissue samples taken from vastly different populations, each with different genetic and environmental backgrounds that may contribute to differing expression profiles. Furthermore, the tissue samples used by each study themselves will be heterogeneous compared to each other. To have adequate quantities of tissue to work with, most studies do high-throughput expression profiling on pooled tissue samples, which results in a gene expression signal that is “averaged” across all cells in the samples (21).

However, the expression of a gene in a single cell may be drastically different from this average. Therefore, depending on how the tissue samples were isolated and which ones were pooled together, the genes may be reported as differentially expressed in an inconsistent direction. One of the studies included in the cancer versus normal comparison (21) investigated the feasibility of pooling tissue samples together by plotting the expression signal of all genes in a pooled sample versus the expression signal of genes from one of the samples in the pool. The authors calculated Pearson correlation coefficients and saw that their values ranged from 0.80 to 0.97, suggesting that the pooling of their specific tissue samples maintained patterns of gene expression representative of each distinct tissue sample. Such an analysis should be done in studies using pooled samples to ensure that the pooled versus unpooled results are comparable.

Finally, poor study design producing inaccurate results may also explain the presence of these genes. In many cases, these genes were ignored because one lone study reported an inconsistent direction of differential expression, which raises suspicions of the validity of the results of the lone disagreeing study. One concern is that some biologically relevant genes may be omitted due to such a study. Therefore, it may be beneficial to include some genes where the majority of the studies agreed on the direction of differential expression instead of the much more stringent “all-or-none” approach we have used. However, because the majority of these genes (75.4%) were reported in only two studies, including these genes would not alter the identity of the highest-ranking candidates greatly (Tables 4 and 5).

Despite of these inconsistencies, we remind the reader that the majority of the multistudy genes (76.6%) were consistently reported as differentially expressed in the same direction, which is an encouraging result, given that each independent study used diverse experimental techniques and tissue samples.

To further assess our results, we performed a literature review of the genes reported by at least seven studies in the cancer versus normal comparison to determine if any have been shown to have diagnostic and/or prognostic utility in colorectal cancer. The most consistently reported differentially expressed gene in our meta-analysis was carbonic anhydrase II (CA2), which was reported as down-regulated in 11 studies. Along with carbonic anhydrase I, CA2 has been shown to have prognostic significance where the expression of both enzymes was related to the metastatic aggressiveness of colorectal cancer (23). Similarly, the potential diagnostic utility of CA2 was shown in a study in which the average level of fecal CA2 in colorectal cancer patients was shown to be significantly greater than those in the control group (24). Immunohistochemistry has been done on colorectal tumor and healthy mucosa tissue to monitor the protein levels of four carbonic anhydrases, among them CA2 (25). That study showed the level of CA2 protein decreased in cancer relative to healthy tissue, thus confirming the transcript based expression profiling results.

Transforming growth factor-β induced 68 kDa (TGFβI) was reported as up-regulated in nine studies. TGFβI is a secreted extracellular matrix protein and was discovered through differential expression analysis of a TGF-β1-treated human lung adenocarcinoma cell line (26, 27). This gene has also been shown to be strongly induced by TGF-β1 in many other human cell lines (28, 29). Despite of the consistent overexpression of this gene, as far as we know, no study has focused specifically on its diagnostic and/or prognostic utility or its role in colorectal cancer tumorigenesis. Overexpression at the protein level has yet to be validated with immunohistochemistry.

IFN-induced transmembrane protein 1 (IFITM1) was also reported as up-regulated in nine studies. IFITM1 has been shown to mediate the antiproliferative properties of the IFN cytokines (30) and was observed to be overexpressed in gastric cancer cells, which resulted in tumor cells being more resistant to natural killer cells and produced a more invasive phenotype (31). As far as we know, immunohistochemistry on human colorectal cancer tissue has not been done for IFITM1 protein; however, reverse transcription-PCR (RT-PCR) was conducted previously on adenomas in a murine model as well as a human colorectal carcinoma cell line, HT29, and elevated expression of IFITM genes (IFITM1, IFITM2, and IFITM3) was observed (32). No further studies have considered the diagnostic and/or prognostic potential of IFITM1 expression in colorectal cancer.

Mal, T-cell differentiation protein-like (MALL), reported as down-regulated in seven studies, is a member of the MAL proteolipid family (33) and encodes an integral protein located in glycolipid- and cholesterol-enriched membranes. To the best of our knowledge, its expression at the protein level has not been measured by immunohistochemistry, and diagnostic and/or prognostic utilities have not been studied.

Carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM1), reported as down-regulated in seven studies, has been shown to be a tumor suppressor in which expression is lost in adenomas and carcinomas. Moreover, the absence of CEACAM1 expression was shown to be correlated with reduced rates of apoptosis in polyps (34). However, a retrospective study did immunohistochemistry on CEACAM1 and showed that 58% of colorectal cancer patients showed an increase in expression (36). It is unclear how the down-regulation of the transcript results in increased CEACAM1 protein expression. Future studies should focus on the half-life of the cancer CEACAM1 transcript to determine if it differs significantly from the normal CEACAM1 transcript. Furthermore, this study did not observe a relationship between CEACAM1 protein levels and overall survival or disease-free survival in colorectal cancer patients (36).

Secreted protein, acidic and rich in cysteine (SPARC) was reported as up-regulated in seven studies and has been shown by our group to be a putative resistance reversal gene (37). Differentially expressed genes between resistant and sensitive human MIP101 colon cancer cells were determined and SPARC was shown to be consistently down-regulated in the resistant cell lines. Their sensitivity was restored by reexpression of SPARC, suggesting that expression of SPARC has prognostic utility. Immunohistochemistry done on colorectal cancer tissue samples showed increased staining of SPARC protein levels (38). However, another immunohistochemistry study (39) showed down-regulation of SPARC due to methylation of its promoter. Further studies related to the role of SPARC in colorectal tumorigenesis are currently underway in our group.

Growth differentiation factor 15 (GDF15), reported as up-regulated in seven studies, is a member of the TGF-β superfamily. Diagnostic and prognostic utility of GDF15 in colorectal cancer has been suggested by studies that showed increased serum levels of GDF15 in colorectal cancer patients relative to healthy controls (40). These levels increased during disease progression and may have clinical use in the management of colorectal cancer patients (41).

MYC, reported as up-regulated in seven studies in our meta-analysis, is a transcription factor that regulates various processes, such as cell cycle progression, differentiation, apoptosis, and cell motility (42). Immunohistochemistry on MYC has shown that its expression increases during disease progression (43), and when combined with nuclear β-catenin expression, MYC expression was shown to have prognostic utility (44).

In conclusion, the results of this meta-analysis identified genes already shown to have diagnostic and/or prognostic potential in colorectal cancer. Perhaps more interesting are the genes, such as TGFβI and IFITM1, which were consistently reported but have yet to be studied specifically as biomarkers. Also, the genes further down the list (that is, those identified as differentially expressed by four, five, six, etc., independent studies) warrant further investigation. Further studies focused on these genes will help in determining a panel of diagnostic and prognostic colorectal cancer biomarkers with sufficient sensitivity and specificity.

The percentage of sequence features that could be mapped to an Entrez Gene ID for each of the three comparison types.

The 77 up-regulated genes reported by three or four studies were further subdivided into semiquantitative categories based on the lowest (Q1), moderate (interquartile range), and greatest (values greater than those in Q3) number of tissue samples to give greater importance to the average fold change criteria for ranking genes when total sample numbers were similar.

The 48 down-regulated genes reported by three or four studies were further subdivided into semiquantitative categories based on the lowest (Q1), moderate (interquartile range), and greatest (values greater than those in Q3) number of tissue samples to give greater importance to the average fold change criteria for ranking genes when total sample numbers were similar.

The 23 up-regulated genes reported by two or three studies were further subdivided into semiquantitative categories based on the lowest (Q1), moderate (interquartile range), and greatest (values greater than those in Q3) number of tissue samples to give greater importance to the average fold change criteria for ranking genes when total sample numbers were similar.

The 16 down-regulated genes reported by two or three studies were further subdivided into semiquantitative categories based on the lowest (Q1), moderate (interquartile range), and greatest (values greater than those in Q3) number of tissue samples to give greater importance to the average fold change criteria for ranking genes when total sample numbers were similar.

Of the 573 multistudy genes, 547 were present in the European Bioinformatics Institute Gene Ontology set of 34,242 annotated genes products. A background set of all genes that were represented at least twice among the platforms used in the cancer versus normal expression profiling studies was used. A total of 24 Gene Ontology terms were found to be statistically overrepresented: 5 biological processes (P), 16 cellular components (C), and 3 molecular functions (F). Number Observed and Total Number show the number of genes from the list found associated with each Gene Ontology term over the total number of genes annotated to that term in Gene Ontology.