Modeling Genomic Diversity and Tumor Dependency in Malignant Melanoma

High-density single-nucleotide polymorphism array hybridization.

Genomic DNA was prepared from near-confluent cells using the DNEasy Tissue Kit (Qiagen) according to the manufacturer’s instructions. For highly pigmented cultures, 1:1 phenol-choloroform extractions were done before column purification. Single-nucleotide polymorphism (SNP) array data were generated using either the StyI chip from the 500K Human Mapping Array set (70 samples; see Supplementary Table S1) or the XbaI chip from the 50K Human Mapping Array set (31 samples; Affymetrix, Inc.). Array experiments using 250 ng of genomic DNA were carried out in a 96-well format using a Biomek FX robot with dual 96 and span-8 heads (Beckman Coulter) and a GeneChip Fluidics Station FS450 (Affymetrix), as described in Supplementary Methods. After hybridization, microarrays were scanned using the GeneChip Scanner 3000 7G (Affymetrix) to generate .cel and .txt files. Subsequent data processing steps are detailed in Supplementary Methods.

High-throughput oncogene mutation profiling and DNA sequencing.

Genomic DNA was genotyped for 238 known mutations in 17 oncogenes as described previously (18). Alternatively, BRAF and NRAS mutations were determined by Sanger sequencing (Agencourt). Available samples were sequenced for phosphatase and tensin homologue (PTEN) coding and splice junction mutations using prevalidated exon primers (Agencourt) and analyzed with Mutation Surveyor software (SoftGenetics, LLC.).

Gene expression profiling.

RNA was purified by a TRIzol extraction protocol (Invitrogen), as described in Supplementary Methods. Affymetrix HT-HGU133A chips were used for gene expression profiling (19). Using 1.0 μg of RNA, synthesis of the cRNA target product, microarray hybridization, and scanning (Affymetrix HT scanner) were done in a 96-well plate format according to the manufacturer’s protocols. Scanned data were processed using MAS 5.0 (Affymetrix) as described in Supplementary Methods.

Bacterial artificial chromosome array data processing.

Bacterial artificial chromosome (BAC) array comparative genomic hybridization (CGH) primary melanoma data (20) were downloaded from the GEO database ({"type":"entrez-geo","attrs":{"text":"GSE2631","term_id":"2631"}}GSE2631). The genomic markers were converted from hg16 to hg17 using the UCSC genome browser.¹¹ Missing data were imputed using the mean of neighboring values. Log₂ ratio values were smoothed using the GLAD segmentation algorithm (21). BAC markers in regions of known copy number polymorphisms were removed¹² before downstream analysis.

Genomic Identification of Somatic Targets in Cancer analysis.

Details of the Genomic Identification of Somatic Targets in Cancer (GISTIC) algorithm are described elsewhere (22). Briefly, each SNP array or CGH data set is defined in terms of copy-number values for N markers along the genome in M tumor samples. To identify significant regions of amplification and deletion, the GISTIC algorithm considers both the frequency F of an amplification or deletion above (or below) a certain threshold defined by normal samples and the average amplitude of that alteration A beyond the threshold. Thereby, a genome significance score (G score) G = F × A is defined over N markers for all M samples. At each individual marker n, the computed metric g is compared with the null hypothesis—a background metric defined by permutation testing across all the markers. Correcting for multiple hypotheses, a false discovery rate (q value) is generated at each marker n, aggregating information from M samples. Using GLAD-segmented copy number data from each sample (21), the minimal common region of overlap is used to define boundaries of peaks in the distribution of q values. A “wide peak” boundary is determined by the samples remaining after leaving out the individual samples that define the borders of the “narrow peak.” This wide peak reduces the influence of individual samples in dictating peak borders. In Tables 1 and ​and22 the wide peak region is indicated.

Inference of loss of heterozygosity.

A Hidden Markov Model (dChipSNP, August 16, 2006 build) that considers haplotype information from a pool of normal samples was used to infer loss-of-heterozygosity (LOH) calls, as previously described (ref. 23; Supplementary Methods).

Hierarchical clustering using genomic alterations.

An input binary matrix was generated from the GISTIC algorithm (see above). For each lesion defined by GISTIC, a ”1” was assigned to a lesion i in a given sample j if the middle marker in the lesion i region exceeded the designated log₂ ratio threshold of 0.3 for amplifications, or −0.3 for deletions. Otherwise, that cell ij was assigned a “0.” Hierarchical clustering (24) with Euclidean distance and complete linkage was applied to this matrix using the GenePattern software package (25).

Significance analysis of microarrays.

Significance analysis of microarray (SAM) was done using TIGR MeV 4.0¹³ to identify differentially expressed genes associated with a broad copy number lesion (i.e., chromosome 10 copy number loss) or a more focal event (i.e., peak of a GISTIC lesion). For the chromosome 10 studies, samples were first ranked into three classes based on mean log₂ copy number over the region. Thresholds were chosen based on the 33rd and 67th percentiles of the copy number distribution. Calculations were done with K-nearest neighbors (10 neighbors), 200 permutations, and the Tusher method (26) to select the minimum δ where the median number of false-positive genes is zero.

Pharmacologic growth inhibition assays.

Cultured melanoma cells were added to 96-well plates at 1,000 (A375 and MCF7) and 3,000 (all short-term melanoma cultures and the MALME 3M cell line) cells per well. Cells were allowed to adhere overnight, and medium containing serial dilutions of CI-1040 (from 100 to 1 × 10⁻⁶ μmol/L; WuXi PharmaTech) was added. Cells were incubated for 96 h, at which time cell viability was measured using the CellTiter-Glo viability assay (Promega). Viability was calculated as a percentage of untreated control values after background subtraction.

“Cross-platform” analyses of melanoma tumor genomic DNA.

Given that the GISTIC algorithm computes genomic regions of significance using smoothed copy number data, we reasoned that this method might also facilitate systematic genomic studies of tumors whose chromosomal copy number data derive from distinct microarray platforms. To explore this possibility, we compared the GISTIC plots generated on 250K SNP arrays from 70 melanoma short-term cultures, 50K SNP arrays from 31 additional short-term cultures and cell lines, 33K tiling arrays from 45 independent cell lines (32), and BAC array CGH data derived from 70 primary cutaneous melanomas (ref. 20; see Materials and Methods; Fig. 1A and ​andB;B; Supplementary Fig. S1). The resulting GISTIC plots were strikingly similar across all data sets, suggesting that this algorithm could extract essential tumor genomic information regardless of the platform used. To confirm these similarities, we carried out pairwise Pearson correlations between GISTIC significance values of cultured melanoma data, primary melanomas, and several additional solid tumor data sets whose chromosomal aberrations had previously been examined (Supplementary Table S2; see Supplementary Methods). Notably, GISTIC correlations between cultured melanoma cells and primary cutaneous melanomas from nonchronic sun-damaged skin were stronger than any other correlation observed among all pairs of primary tumor sets, approaching that of the positive control (see Supplementary Methods). Thus, the chromosomal regions of significance gleaned from GISTIC analyses were robust to the types of microarray platform variances common to cancer genomic data sets.

The genomic diversity of melanoma in vitro and in vivo.

Theoretically, cultured human cancer cell collections may represent a biased malignant subset that fail to encompass the pathophysiology relevant in primary tumors. To address this at a genomic level in cutaneous melanoma, we compared the genomic regions of significance found in cultured melanoma cells to those observed in primary tumors. As shown in Fig. 1A and ​andB,B, the landscape of genomic alterations was strikingly similar between cultured cells and primary tumors, even with respect to subtle contours within the GISTIC plots. Some significant alterations present in cultured cells were absent in primary tumors [e.g., gains on chromosome 3p14 (MITF) and chromosome 15q, and a region of loss on 4q]; in most cases, these were attributable to low BAC clone coverage in the primary tumor data. Conversely, a single peak of chromosomal deletion centered at 13q33 exhibited reduced prominence in cultured cells compared with primary tumors (Fig. 1B). Interestingly, the underlying 3.2-Mb locus spanned the ERCC5 excision repair gene, and this deletion showed increased significance in primary melanomas from chronic sun-damaged skin (ref. 20; Supplementary Fig. S2). These findings raised the possibility that deletions involving ERCC5 may contribute to melanoma genesis in the setting of chronic sun damage. Overall, this comparative analysis suggested that, in aggregate, cultured melanoma collections might provide a genetically appropriate model system for systematic functional genomic characterization.

Integration of significant genomic lesions with melanoma gene expression data.

Although the majority of GISTIC peaks pinpointed genomic regions implicated previously in melanoma (refs. 20, 32, 33; Supplementary Fig. S1A), in most cases the relevant effector genes remain unknown. To identify candidate genes targeted by the GISTIC lesions, we systematically partitioned our data set into two classes based on the presence or absence of the relevant GISTIC lesion and used the resulting classes to carry out SAM analysis of sample-matched gene expression data (see Materials and Methods). In principle, amplified genes that exhibit significantly up-regulated expression may denote candidate oncogene targets, whereas those within regions of deletion that show significant down-regulation may exert tumor suppressor mechanisms.

The results of this genome-wide SAM analysis are shown in Tables 1 and ​and2.2. In several instances, a short list of differentially expressed genes was identified whose cellular functions suggested plausible cancer mechanisms. For example, we identified five up-regulated genes from the 15q26 GISTIC region that included IQGAP1 (Table 1), which encodes a GTPase-activating protein that modulates RAC1/cdc42-mediated cellular adhesion (34). Similarly, four up-regulated genes from the GISTIC locus at 22q13 included the EP300 gene (Table 1), which encodes a key transcriptional coactivator that cooperates with numerous cancer-associated factors, including microphthalmia-associated transcription factor (MITF; ref. 35). Several known or candidate melanoma tumor suppressors were also suggested by this analysis (Table 2), including PTEN (10q23) and KLF6 (10p15; described below). Thus, combined analyses of sample-matched, genome-scale data posit multiple candidate tumorigenic mechanisms enacted by melanoma genomic aberrations.

Genome-wide measurements of LOH in melanoma.

High-density SNP arrays also enable inference of LOH (36), even in the absence of matched normal DNA (23). To identify significant LOH regions within our cultured melanoma collection, we predicted LOH patterns within each sample using a hidden Markov model (23), and applied the GISTIC algorithm to the resulting LOH segments (see Materials and Methods). A strong concordance was observed between regions of LOH and copy loss (Fig. 1C), suggesting that most melanoma LOH events occur by hemizygous deletion. The major exceptions to this pattern occurred at chromosomes 5q and 17p, where LOH in the absence of significant hemizygous deletion (“copy-neutral” LOH) was apparent. Although the relevant tumor suppressor at chromosome 5q has not been identified, the LOH peak at chromosome 17p centers on the well-known TP53 tumor suppressor gene (Fig. 1C). This suggests that copy-neutral LOH may represent an important mechanism of p53 inactivation in melanoma. In a subset of short-term culture samples, a segment of chromosome 4q contained hemizygous deletion but limited LOH (Fig. 1C). This unusual phenomenon, which could potentially be explained by copy loss and biallelic retention in the setting of hyperploid tumors, implies that this region of the melanoma genome might be relatively refractory to LOH or uniparental disomy.

Molecular classification of melanoma based on chromosomal aberrations.

We next sought to determine if significant chromosomal events might organize the melanoma samples into molecularly distinct and/or biologically relevant subsets. To address this, we applied established unsupervised learning methods (37, 38) to the lesions identified by GISTIC (see Materials and Methods). As shown in Fig. 2, hierarchical clustering grouped GISTIC-annotated cultured melanoma samples into six main branches. Two major subclusters (clusters 4 and 5) were characterized by the presence of chromosome 7p/7q gains (P = 1.1e–05, P = 2.5e–05, two-tailed Fisher’s exact test) and chromosome 10 losses (P = 7.9e–4, P = 6.5e–05) with cluster 5 being enriched for additional alterations, most notably gains on chromosome 22q (P = 1.4e-3). These clusters were associated with an enrichment of BRAF mutation compared with NRAS mutation (P = 4.7e-3, two-tailed Fisher’s exact test; Fig. 2), thereby highlighting the observation that in a subset of samples mutated BRAF is amplified, at least to low levels, in melanoma (Supplementary Fig. S3C–E). Another prominent subcluster (cluster 6) was largely defined by a relative paucity of chromosome 10 losses when compared with clusters 1 to 5. Interestingly, although 7q gains were relatively common in this terminal branch, the prevalence of BRAF(V600E) mutation was diminished compared with clusters 4 and 5, raising the possibility of an additional target(s) of 7q copy gains in melanoma. In contrast, clusters 1 and 2 exhibited a reciprocal pattern, with chromosome 10 losses characteristic of both branches. Overall, these results suggested that clustering GISTIC lesions may provide biologically relevant molecular groupings of melanoma tumor samples.

Open in a separate window

Figure 2.

Clustering analysis of significant melanoma genomic alterations. Hierarchical clustering of GISTIC lesions (discretized smoothed copy number; see Materials and Methods) by the Euclidean distance metric and complete linkage. Rows, genomic lesions identified by GISTIC algorithm (see text); columns, samples. Red, presence; blue, absence of lesions denoted by cytobands (A, amplification; D, deletion; right). Mutation status of BRAF and NRAS are noted above the matrix. Major clusters are indicated by boxes labeled 1 to 6.

Predicting genes targeted by nonfocal chromosome 10 deletions in melanoma.

In melanoma, as in many solid tumors, most statistically significant chromosomal aberrations span large genomic regions whose target genes are unknown (20, 32, 33). Among many other lesions, this pattern is exemplified by the gains and losses involving chromosomes 7 and 10, respectively, two prominent drivers of the GISTIC-based melanoma subclusters described above. Although one of the most significantly “amplified” regions on chromosome 7 centers at 7q34 (contains BRAF; Fig. 1A), most chromosome 10 copy losses cover the entire chromosome, punctuated by a more focal GISTIC peak at 10q23.31 (Fig. 1B; Supplementary Fig. S1A and B) spanning the PTEN tumor suppressor gene. These initial observations are consistent with a model in which chromosome 7 gains and chromosome 10 deletions enhance oncogenic BRAF and diminish PTEN expression, thereby suggesting cooperating effects of these two cancer genes in melanoma, as previously suggested (39). On the other hand, whereas broad chromosome 10 deletions are highly prevalent even in early-stage primary melanomas (ref. 20; Fig. 1B), the relative contribution of PTEN inactivation in melanoma establishment remains unclear.

To investigate this, we examined the extent to which panchromosome 10 deletions were associated with functional PTEN inactivation, either through genetic mutation or mRNA/protein loss. As shown in Fig. 3A, 5 of 42 short-term cultures and cell lines harboring large chromosome 10 losses also contained focal PTEN homozygous deletions; similar results were observed in an independent melanoma cell line panel previously described (Supplementary Fig. S3A). However, consistent with prior DNA sequencing analyses of short-term melanoma cultures and primary tumors (40, 41), only 8 of 78 lines examined by sequencing contained inactivating PTEN mutations; of these samples, only 5 also harbored hemizygous loss and LOH of chromosome 10. Using sample-matched gene expression data, we then examined PTEN mRNA expression in relation to genomic deletion. Here, whereas PTEN homozygous deletion correlated with a marked reduction in PTEN expression, <50% of samples with hemizygous deletions exhibited decreased PTEN levels when compared with samples without alterations at this locus (Fig. 3B). Similar results were observed following immunoblot analysis of PTEN protein levels in a subset of samples (Fig. 3C). Although tempered by the incomplete sensitivity to detect small-scale deletions at the PTEN locus, these observations raised the possibility that an as yet uncharacterized tumor suppressor gene(s) located on chromosome 10 might contribute to melanoma establishment in addition to PTEN.

Open in a separate window

Figure 3.

Integrative analysis of chromosome 10 deletions in melanoma. A, heatmap view of smoothed 250K and 50K SNP array data spanning the PTEN locus (red, copy gains; blue, copy losses). Samples (columns) are sorted based on copy number values derived from segmented data. B, PTEN gene expression values stratified according to homozygous deletion, hemizygous loss, or retention of the underlying locus. C, summary of PTEN protein levels in relation to chromosomal deletions affecting the PTEN locus. D, significance analysis of melanoma gene expression data. Samples are stratified according to loss or retention at chromosome 10 (see Materials and Methods). The positions of CUL2, KLF6, and PTEN are shown.

To identify candidate tumor suppressor genes enacted by nonfocal chromosome 10 deletions in an unbiased manner, we partitioned 66 matched samples into tertiles based on the ranked magnitude of chromosome 10 loss (see Materials and Methods) and segregated the two extreme tertiles according to the class distinction “chromosome 10-deleted” (22 samples) versus “chromosome 10-wild-type” (22 samples). We then carried out a genome-wide significance analysis (SAM) using sample-matched gene expression data, as described above. As shown in Fig. 3D, 43 unique transcripts exhibited significantly diminished expression in association with chromosome 10 deletion at a δ value of 1.4673 (median false number of genes = 0). Interestingly, two chromosome 10 genes with putative tumor suppressor roles—CUL2 (chromosome 10p11) and KLF6 (chromosome 10p15)—scored even more highly than PTEN in this analysis (Fig. 3D; Supplementary Fig. S6). CUL2 encodes a cullin protein whose yeast homologues negatively regulate the cell cycle (42) and also forms an E3 ubiquitin ligase complex with the von Hippel-Lindau tumor suppressor protein (43, 44). KLF6 encodes a zinc finger transcription factor implicated as a tumor suppressor in several malignancies (45–47). Notably, the KLF6 locus was also identified as significantly deleted by GISTIC analysis (Table 2). Thus, CUL2 and KLF6 represent intriguing candidate melanoma tumor suppressor genes targeted by chromosome 10 deletions.

The Swiss National Science Foundation grant 32-63704.00; Gottfried and Julia Bangerter Rhyner Stiftung; Genentech, Inc.; National Cancer Institute grants K08CA115927 and P50CA93683; the Burroughs-Wellcome Fund; and the Robert Wood Johnson Foundation. T.R. Golub is an investigator of the Howard Hughes Medical Institute.