Gene Expression Changes in an Animal Melanoma Model Correlate with Aggressiveness of Human Melanoma Metastases

Introduction

Metastasis is the dispersal of cancer cells from their primary loci to distant organs, and accounts for more than 90% of deaths in cancer patients. The mechanisms of metastasis remain incompletely understood (1-3). Metastasis is a rare event, as shown by both clinical and animal studies. Clinicians often find a large number of circulating tumor cells in cancer patients, but not many detectable metastases (4, 5). The process of metastasis has also been studied in animal metastasis models. In one such model, a pool of poorly metastatic human tumor cells is injected into immunodeficient mice and the resulting metastases are isolated and cultured in vitro as cell lines (6-10). These cell lines often show enhanced metastatic ability when reinjected into immunodeficient mice. Studies using such models have found that most of the injected cancer cells are able to disseminate into different organs but only a small proportion of them grow as detectable metastases, consistent with the notion that metastasis is a relatively rare event. Consequently, metastases have been postulated to result from small populations of cancer cells within a primary tumor, which are able to enter and survive in the circulation and then exit the circulation and grow in a distant organ (4, 5, 11).

Recently, microarray analyses have provided valuable tools to dissect further the mechanisms of cancer progression and to improve cancer treatment. Gene expression patterns have been used to classify the subtypes of primary tumors and to predict the clinical outcome of their treatments (12-14). Among these studies, several have shown that primary tumors and their metastases have similarities in gene expression profiles (15-17), and van't Veer et al. (13) and Ramaswamy et al. (14) discovered that some primary tumors contain gene signatures that can predict their propensity to metastasize. Those observations raise the possibility that the ability to metastasize is determined early in primary tumor development and does not require further selections among primary tumor cells. They suggest an alternative to the view that metastases arise from rare populations within the primary tumor.

These two views of metastasis are not mutually exclusive (18). It is entirely possible that different primary tumors can be, either ab initio or as a consequence of their progression/evolution, of good or poor prognosis and for these properties to be reflected in their gene expression profiles. That does not exclude the possibility that further alterations in gene expression can be either contributory to or necessary for effective spread and growth of metastases, and indeed, data are available showing the existence of gene expression signatures characteristic of metastases (7, 9, 14). The challenge is to determine gene expression signatures that contribute to various aspects of tumor progression (predisposition to metastasis, actual metastasis, aggressiveness of metastases, etc.) and to relate those signatures to clinical data and outcomes.

We describe here the derivation of a series of metastatic human melanoma cell lines from poorly metastatic parental lines and their analysis in xenotransplant tumor and metastasis models. Genes differentially expressed between tumor samples derived from highly metastatic derivatives and from their poorly metastatic parental lines were identified. Expression of this gene expression signature in human metastases was found to correlate with poor survival of melanoma patients with metastases. Among these genes, many are secreted or membrane proteins, suggesting the importance of interactions between tumor cells and their microenvironment in the aggression of metastases.

Results

Derivation of Two Groups of Highly Metastatic Cell Lines from Poorly Metastatic Human Melanoma Cells

To study the mechanisms of metastasis, we took advantage of an experimental metastasis assay to derive several melanoma cell lines from two poorly metastatic parental lines (Fig. 1A and B ; ref. 8). The parental lines were either from American Type Culture Collection or from Dr. Isaiah Fidler (The University of Texas M. D. Anderson Cancer Center, Houston, TX). These are both A375 melanoma cells but have been cultured separately for decades. Nevertheless, they both remain poorly metastatic when tested in vivo and give rise to few lung metastases when injected into the circulation of immunodeficient mice. These few lung metastases were isolated and amplified in vitro as cell lines (MA-1, MB-1, and MC-1 and MEA-1, MEC-1, and MED-1; Fig. 1A and B). The derivatives from the ATCC parental line were denoted as set A cells, and those from the Fidler parental line were denoted as set F cells. These cells were reinjected into mice for a second round of selection to generate MA-2 and MC-2 from set A and MEA-2, MEC-2, and MED-2 from set F (Fig. 1A and B). The MEA-2 cell line was injected once more into immunodeficient mice to generate cell line MEA-3 (Fig. 1B). The derived cell lines exhibit increased metastatic ability as compared with their parental lines when tested by i.v. injection (Fig. 1C).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 1.

Derivation of highly metastatic melanoma cell lines. A and B. Lineages of derived cell lines. Two related poorly metastatic A375P melanoma cell lines were injected into immunodeficient mice and individual lung metastases were isolated and cultured in vitro as independent cell lines (see Materials and Methods). C. These derived cell lines are more metastatic than the parental lines. For set A cell lines, 2 × 10⁵ cells were injected into immunodeficient mice (n > 5) and lung metastases were counted after 1 mo. For set F cell lines, 5 × 10⁵ cells were injected and lung metastases were counted after 2 mo. *, P < 0.05.

Tumors from the Two Groups of Metastatic Derivatives Show Similar Expression Profiles

We did microarray analyses to identify genes that are up-regulated or down-regulated in the tumor samples from highly metastatic cells compared with those from the parental lines. Each cell line was injected into immunodeficient mice either intravenously or subcutaneously. The lung metastases or subcutaneous tumor from each mouse were processed as one biologically independent sample in the microarray analyses after their RNA was extracted and hybridized onto human oligonucleotide microarrays (see Materials and Methods). Genes differentially expressed between the tumor samples from metastatic variants and those from their respective parental lines were identified using GenePattern software. The expression levels of some of them (∼30 genes) were validated by real-time PCR and found to be largely consistent with the results from array analyses. The top and bottom 50 differentially expressed genes are shown in Fig. 2 . Gene set enrichment analysis (GSEA) was done to examine the similarities between the expression profiles (19). GSEA uses an algorithm that measures the cumulative enrichment of one set of genes (called a gene set) in a ranked second gene list from an array comparison. In our analysis, the most up-regulated or down-regulated 100, 200, or 500 genes from the tumor samples of one set of metastatic cells were used as a gene set to measure their enrichments in the tumor samples from the other set. A normalized enrichment score was assigned to each gene set and the statistical significance of its enrichment was measured by P values and false discovery rate (FDR) scores (see Materials and Methods). The results showed that genes altered in the tumor samples from one set of cells were significantly altered in the other (P < 0.05 and FDR < 0.25; Table 1 ), indicating that similar genes or pathways are regulated in the tumor samples from the two independent sets of derived human metastatic melanoma cells, although the precise order of genes differed between the two sets. However, although the gene lists from the two sets of tumors were related, the gene expression values for set A samples were more consistent and the set A data set has been used in the subsequent analyses.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2.

GSEA shows that data sets A and F contain similar patterns of gene expression. Expression patterns of the top and bottom 50 differentially expressed genes in set A and set F data identified by conventional marker selection method (class neighborhood in GenePattern software). Red boxes, samples from the poorly metastatic parental lines. Sample names and gene symbols on each heat map are also listed in Supplementary Document S2 for reference.

Metastasis Genes from Set A Tumors Are Differentially Expressed in Human Melanoma Metastases and Associated with Poor Survival

We next tested whether the up-regulated or down-regulated genes in tumor samples from our metastatic variants correlate with the clinical outcome of human melanoma patients. We selected probe sets that are up-regulated or down-regulated in the tumors from set A metastatic cells by at least 3-fold and with a maxT value <0.01 (see Materials and Methods). One hundred eighty-five probe sets (99 up-regulated and 86 down-regulated) were selected, which correspond with 150 nonredundant genes (74 up-regulated and 76 down-regulated; Supplementary Document S3). These probe sets were used as templates for the following nearest template prediction method (for details, see Materials and Methods). Briefly, the expression values of these probe sets in array data from 52 human melanoma metastases were extracted and normalized with a mean equal to 0 and a SD equal to 1. Their distances from the templates were calculated using Pearson's correlation coefficient and FDR scores were calculated to measure the statistical significance. The human metastasis samples were separated into three major classes based on their distances from the templates and their FDRs (Fig. 3A ). Class 1 and class 3 metastasis samples show opposite expression patterns of the 185 probe sets and their FDR scores are <0.05. Class 1 samples showed similar expression patterns to the tumor samples from our set A metastatic derivatives, whereas class 3 samples show similar expression patterns to those from our set A parental line. Class 2 samples have FDR scores of >0.05 and their expression pattern of the 185 probe sets is in between class 1 and class 3 samples. The survival data from patients with the three classes of metastases were used to generate Kaplan-Meier survival curves (Fig. 3B; Supplementary Document S4). Patients with class 1 metastases had shorter survival than patients with class 2 or class 3 metastases. The difference in survival probability between patients with class 1 and class 3 metastases was found to be significant by a log-rank test (P = 0.003). Because class 3 metastases include a few lymph node metastases, which contain high percentages of contaminating lymphocytes, we tested whether the survival difference between class 1 patients and class 3 patients still exist when the lymph node metastases were excluded from our analyses. The results showed that the survival probability of class 1 patients still differed significantly from that of class 3 patients when lymph node metastases were omitted from the analyses (Supplementary Fig. S5). These data suggest that the set of 150 genes differentially regulated in the tumors from our set A metastatic cells correlates with the aggressiveness of human melanoma metastases.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3.

Genes up-regulated in the metastatic derivatives from set A correlate with metastatic death in human melanoma patients. A. The 185 probe sets (150 genes; see text) were used as the template for the nearest template prediction (for details, see Materials and Methods). They include 99 up-regulated probe sets (orange bar) and 86 down-regulated probe sets (green bar). Expression values of the 185 probe sets were extracted from expression profiling data from human melanoma metastases and were used to calculate the distance of each sample from the template. The human melanoma metastases were separated into three classes based on their distance and FDR value. Class 1 metastases (red bar) express high levels of the up-regulated genes in the signature, FDR < 0.05; class 3 metastases (blue bar) express high levels of the down-regulated genes in the signature, FDR < 0.05; and class 2 metastases (gray bar) express intermediate levels of the signature genes, FDR > 0.05. B. Kaplan-Meier survival curves were generated based on the correlation between the survival of patients and the classes of metastases they carry. The difference in survival probability between patients with the class 1 and class 3 metastases was found to be significant by a log-rank test. C. Similar analyses were done as in A, but using expression data from primary human melanomas. The 185 probe sets (orange bar, up-regulated probe sets; green bar, down-regulated probe sets) separated the primary melanomas into three classes (red, gray, and blue bars, respectively). Black dots, those melanomas that gave rise to metastases. Their Breslow thicknesses and Clark levels are labeled by colored hexagons and squares, respectively. The two extreme classes (classes 1 and 3) of primary melanomas do not differ in their abilities to develop metastases (P = 0.224, Fisher's exact test), their Breslow thicknesses (P = 0.32 if the proportion of samples with thickness of <1 mm was compared between the two classes, Fisher's exact test), or their Clark levels (P = 0.38 if the proportion of level III samples was compared between the two classes, Fisher's exact test). They also did not differ significantly in survival probability (P = 0.106, log-rank test), possibly because of the small number of patients available for analysis (see Discussion).

Interestingly, expression levels of these 150 genes also separated 31 human primary melanomas into three classes (Fig. 3C). However, because the nearest template prediction method describes the relative distribution among the input samples, this classification could be specific to primary melanomas. To test whether the three classes of primary melanomas correspond with the three classes of metastases as described in Fig. 3A and C, we mixed the primary melanoma samples with metastasis samples and regrouped them using the nearest template prediction method. We found that most of the class 1 primary melanomas grouped with class 1 metastases, and most of the class 3 primary melanomas grouped with class 3 metastases (Fig. 4, bottom bars ), suggesting that the gene expression profile characteristic of aggressive human melanoma metastases (see Fig. 3A and B) is already represented in a subclass of primary melanomas (Fig. 3C). However, patients with class 1 and class 3 primary melanomas did not differ significantly in their survival (P = 0.182, log-rank test; Supplementary Fig. S6). They also did not differ in their Clark levels when the proportion of level III samples in all class 1 samples was compared with that in all class 3 samples (Fig. 3C, colored squares; P = 0.38, Fisher's exact test) or in their Breslow thicknesses when the proportion of samples with thickness of <1 mm in all class 1 samples was compared with that in all class 3 samples (Fig. 3C, colored hexagons; P = 0.32, Fisher's exact test). They also did not differ in their propensity to develop metastases (Fig. 3C, black dots; P = 0.224, Fisher's exact test). This lack of significant difference between these two classes of primary melanomas may well be due to the small number of patients that are available for the analyses (see Discussion).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4.

The three classes of primary melanomas and metastases show similar expression patterns of the 150 genes. The expression values of the 150 genes (see text) were extracted from the mixed primary melanoma and metastasis samples and these samples were separated into three classes by the nearest template prediction method as described in Fig. 3A. These classes were indicated by red, gray, and blue bars on the top of the heat map. Each class of primary melanomas and metastasis samples was assigned with the same class number as in Fig. 3, indicated by the red, gray, and blue bars at the bottom of the heat map.

The correlation between our 150 genes and the aggressiveness of human melanoma metastases was also tested by GSEA. We first preranked the ∼22,000 probe sets on the HU133A DNA chip based on their Cox-ranking scores, measuring their correlations with the survival of patients carrying melanoma metastases. Then the 185 probe sets that correspond with our 150 genes were split into two gene sets: up-regulated genes in one set and down-regulated genes in the other. The enrichments of these two gene sets in the above preranked probe set list were measured by GSEA. The results showed that the up-regulated probe sets were significantly enriched (P < 0.05 and FDR < 0.25; Table 2 ; Fig. 5 ; Supplementary Document S7), and 56 of them contributed positively to the enrichment (Supplementary Document S7). These data suggest that up-regulated genes among the 150 correlate with poor survival of patients carrying melanoma metastases. This is consistent with our previous finding using the nearest template prediction method and suggests that the up-regulated genes are most significantly correlated with aggressiveness in melanoma metastases.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5.

GSEA shows that up-regulated genes of the 150 correlate significantly with poor survival of patients carrying metastases. The ∼22,000 probe sets on HU133A Affymetrix DNA chips were ranked by Cox scores based on their correlation with the survival of the melanoma metastasis–bearing patients. This preranked list of probe sets was used as the template to assess the enrichments of the 99 up-regulated probe sets and 86 down-regulated probe sets from our signature. The up-regulated probe sets (74 genes) were found to be significantly enriched and to correlate with poor survival of the metastasis-positive patients. The enrichment of the up-regulated gene set is shown schematically. The X-axis of the curve for enrichment scores includes the ∼22,000 probe sets on the HU133A chip, with those correlating best with poor survival on the left and those correlating best with good survival on the right (shown in the Ranked List Metric below the curve). Each probe in the gene set is shown as a vertical line underneath the X-axis of the curve, and the cumulative enrichment score is plotted as the green curve reaching a maximum enrichment at a score of 0.54. The rank order and the contribution of each of the 74 up-regulated genes to enrichment (shown as Core Enrichment) are listed in Supplementary Document S7.

Strikingly, among these 150 genes, many are externally exposed proteins, such as secreted or extracellular matrix proteins or plasma membrane proteins (Supplementary Document S1 and Fig. 6B ): these two categories represent about 56% of the 74 up-regulated genes and 42% of the 76 down-regulated genes compared with ∼15% in the total genes unselected by fold changes and maxT value (which correspond to 9,108 probe sets on HU133A chip; Fig. 6A; Supplementary Document S8). KEGG pathway analyses using DAVID/EASE software showed that the 150 genes seem to be involved mostly in processes such as extracellular matrix-receptor interaction, focal adhesion, transforming growth factor-β signaling pathway, as well as glycerolipid metabolism and arginine and proline metabolism (Supplementary Document S8). For further discussion on the potential implications for metastasis of this gene signature, see Discussion.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6.

The 150 genes include many genes encoding secreted proteins. A. Before applying the selection criteria of a fold change >3 and maxT <0.01, all the genes that passed the preprocessing step (corresponding to 9108 probe sets in total; see Materials and Methods) were categorized based on their known or predicted cellular distributions in the GeneOntology database. B. The 150 genes were categorized based on their known or predicted cellular distributions in GeneOntology database. External proteins, including secreted extracellular and membrane proteins, represent 56% of the up-regulated genes and 42% of the down-regulated genes.

Discussion

In this article, we report the derivation of highly metastatic human melanoma cell lines from poorly metastatic parental lines using an animal metastasis model. We subsequently identified a “metastasis aggressiveness gene signature” by comparing the gene expression patterns of tumor samples from the highly metastatic derivatives with those from their parental lines. By comparisons with gene expression data from human clinical samples, we found that expression of this metastasis gene signature in human melanoma metastases correlates with poor survival of the corresponding patients. The signature is able to segregate melanoma-bearing patients into three groups, one of which has a significantly lower survival probability. This suggests that the signature provides an indication of “aggressiveness” of the melanoma metastases rather than of metastasis per se, similar to the lung metastasis signature reported by Minn et al. (10). This result has been confirmed by alternative methods such as GSEA and hierarchical clustering (data not shown). Interestingly, our gene signature is also able to separate primary tumors into the same three classes, a result reminiscent of some other gene signatures that have been reported (13, 14, 17). Given the size of the sample of primary tumors, we were unable, as yet, to show a predictive role for the gene signature when detected in primary tumors. Only six patients' primary tumors expressed the signature (see Fig. 3C, class 1), and in those six patients, only one had reported metastases and two died—insufficient for statistical analysis. It will be of interest to test this gene signature on independent data sets from larger numbers of patients to test whether it has useful predictive value. It is clear that not all melanoma metastases express the 150-gene signature that we have described, and the data suggest that those that do express it have a poorer prognosis, so the signature may well have prognostic value for patients with diagnosed and biopsied metastases. Those possibilities need to be tested against independent sets of melanoma samples.

Turning next to the implications of our results for understanding the cellular mechanisms of metastatic spread, several issues need to be discussed. Because the 150-gene signature is expressed in some primary tumors, it falls into the category of expression profiles that preexist in the bulk of (some) primary tumors, although it is unclear whether the profile that originates in the cells of origin of the tumor is a consequence of the nature of the initial oncogenic transformation event or arose during progression/development of the primary tumor. It is also unclear whether this signature is necessary or sufficient for metastasis, although the data do suggest that it correlates with, and may contribute to, the aggressiveness of the metastases. It is entirely possible that other sets of genes are necessary for metastasis even in the context of cells expressing the 150-gene signature, and it is clear (from the existence of metastases that do not express this signature) that it is not the only gene expression state characteristic of metastasis.

It is of some interest that the 150-gene signature that we report here is significantly enriched in secreted and membrane proteins that could be involved in tumor-microenvironment interactions contributing to the progression of metastases. Among them, the 56 genes that contribute positively to the enrichment of the 74 up-regulated genes in aggressive melanoma metastases are likely good candidates for further investigations. Many of them have previously been implicated in melanoma progression—examples include endothelin receptor (20), ERBB (21), Frizzled homologue 7 (22). All these genes can be tested directly in the animal metastasis model for their roles in the interactions between tumor cells and their microenvironment in metastasis. They should also be good markers for diagnosis and prognosis of malignant melanoma because they are secreted and thus may be present in the plasma of cancer patients and easily detected by methods such as ELISA or mass spectrometry (23).

In conclusion, the findings we report here should contribute to the understanding, diagnosis, and treatment of malignant melanoma, which currently remains essentially untreatable.

Materials and Methods

Derivation of Metastatic Melanoma Cell lines

The highly metastatic human melanoma cell lines were derived from two poorly metastatic parental lines as described (8). The two parental A375 lines were obtained either from American Type Culture Collection (ATCC CRL-1619; for set A cells) or as a gift from Dr. Isaiah Fidler (for set F cells). Briefly, 500,000 A375 cells were injected i.v. into nude mice (Cby.Cg-Foxn1^nu, The Jackson Laboratory). Two months later, individual lung metastases (believed to be clones) from different mice were harvested and amplified in vitro as independent cell lines. These cell lines were reinjected into mice for a second round of selection. MEA2 cell line was reinjected for a third round of selection. In total, four cell lines were derived from the set A parental line and seven were derived from the set F parental line (Fig. 1A and B). An SM cell line was derived in a similar manner in Dr. Fidler's laboratory (24) and was also included in the array analyses along with the other set F cell lines.

To test the metastatic properties of derived cell lines, 200,000 set A cells or 500,000 set F cells were injected i.v. into nude mice and lung metastases were counted under a dissecting microscope 2 mo later. These lung metastases were also harvested for microarray analyses. Those from the same mouse were pooled and processed as one sample for microarray analyses. To generate subcutaneous tumors, 500,000 cells were injected subcutaneously into the right flank of immunodeficient mice and the tumors were harvested 4 wk after injection. Each subcutaneous tumor was processed separately and used as one sample for microarray analyses. All the samples included in the array analyses and their corresponding cell lines were listed in Supplementary Document S1.

RNA Preparation and Data Collection

Each cell line was injected into at least three different mice to obtain subcutaneous tumors or lung metastases for microarray analyses. As described above, tumors from the same mouse were pooled and processed as one sample for the array analyses. RNA was extracted from the tumors using Qiagen RNeasy Midi Kit according to the manufacturers' instructions. cRNA was prepared according to the GeneChip Expression Analysis Technical Manual (Affymetrix), hybridized onto HU133A chips (Affymetrix), and scanned by a GeneArray@ 2500 Scanner (Affymetrix).

The quality of raw microarray profiles was generally assessed by measurements of overall microarray fluorescence intensity (e.g., mean, variance), the distribution of feature or spot intensities, and the proportion of total genes showing significant signal. Thirty-two set A and 39 set F data sets passed these quality control criteria and were used for subsequent data analyses such as normalization, expression marker selection, and GSEA. The data have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (GEO)⁷ and are accessible through GEO Series accession nos. GSE7929 (set A) and GSE7956 (set F). They were normalized using RMAexpress software⁸ before being further analyzed.

Eighty-three fresh melanoma biopsies from patients undergoing surgery were collected from 1992 to 2001 as a part of the diagnostic workup or therapeutic strategy. Immediately after surgery, half of each specimen was fixed in formalin and processed for routine histology, and the other half was immediately snap-frozen and stored in liquid nitrogen until use for RNA extraction. Histopathologic diagnosis of each tissue specimen was done independently by two histopathologists. All patient material has been collected and used according to the approval by the institutional ethics committee and written informed consent obtained in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki. Total cellular RNA was prepared by guanidinium thiocyanate extraction and cesium chloride centrifugation and purified from remaining melanin with the Qiagen RNeasy Fibrous Tissue Mini Kit. cRNA was prepared and profiled as described above. The data have been deposited in the National Center for Biotechnology Information GEO⁷ and are accessible through GEO Series accession no. GSE8401.

Microarray Analyses

Marker selections were done using the GenePattern software⁹ (25). Set A data were first preprocessed using the PreProcessDataset module (after which 9,108 probe sets are left for further analyses) and then marker genes were selected using the ComparativeMarkerSelection module. Genes differentially regulated in tumors from set A metastatic cells were selected based on their fold change (>3) and their adjusted P value (max T < 0.01). They include 185 probe sets, which correspond with 150 nonredundant genes. The GeneOntology analyses for their cellular distribution and the pathways they are involved in were done using DAVID/EASE software.¹⁰

S-plus software (Insightful) was used to generate Kaplan-Meier survival curves and perform the log-rank tests for assessment of statistical significance between each pair of curves. All the Fisher exact tests were done using the web-based calculation tools.¹¹

GSEA

For the analyses shown in Fig. 2, the 100, 200, or 500 most up-regulated or down-regulated genes in tumors from each set of metastatic derivatives were selected using the ClassNeighborhood module in the GenePattern software. Signal-to-noise ratio was used to calculate statistics and the cutoff P values were assigned based on 1,000 random permutation tests. The selected marker genes from one data set were subsequently used as gene sets to measure their enrichment in the other data set by GSEA¹² (see ref. 19). A normalized enrichment score was calculated based on the size of the gene set and its enrichment score. A nominal P value was calculated after permutation testing of the microarray samples and a FDR (26) was calculated to correct for multiple hypothesis testing. Generally, a gene set is considered significantly enriched when its P value is <0.05 and FDR score is <0.25 (19).

For the analysis shown in Fig. 5, the ∼22,000 probe sets on HU133A Affymetrix DNA chips were ranked by their Cox scores, which evaluate the correlation of their expression values with the poor survival of patients that developed melanoma metastases. This preranked list of probe sets was used as the template and up-regulated or down-regulated genes within the 150-gene signature were used as gene sets to perform GSEA and measure their correlation with poor survival of patients with metastases.

Nearest Template Prediction Method

The nearest template prediction is a variation of the k-means clustering. First, we defined the “templates” of the “metastasis” and “nonmetastasis” patterns, which are equivalent to the centroid in the k-means method. They include the 185 probe sets of our signature. In the templates, only the information on the direction of gene expression change was retained [i.e., in the metastasis template, values for the up-regulated probe sets (99 in total) were set to 1, and the values for the down-regulated probe sets (86 in total) were set to 0, and vice versa for the nonmetastasis template].

The expression values of each of the 185 probe sets were normalized across all the samples, with a mean equal to 0 and SD equal to 1. The distance of each sample to either of the templates was assessed by Pearson correlation coefficient. The samples were then separated into three major groups: first based on which template they were closer to, then by the significance of their proximity. The significance of the proximity was evaluated based on an empirical null distribution of the correlation coefficient generated by randomly picking the same number of genes from the entire microarray data for each sample (n = 1000). A nominal P value was computed using the rank of the observed correlation coefficient in the null distribution. The nominal P value was corrected for multiple hypotheses testing using the FDR. An FDR <0.05 was regarded as significant. Samples closer to the metastasis template and with FDR <0.05 were grouped as class 1; samples closer to the nonmetastasis template and with FDR <0.05 were grouped as class 3; and the rest of samples were grouped as class 2.