A new 12-gene diagnostic biomarker signature of melanoma revealed by integrated microarray analysis

Microarray datasets

This study examines the differential expression of genes between normal skin and/or benign nevi, and metastatic melanoma using a meta-analysis approach. The experimental protocol of this study is shown in Fig. 5 and commenced with the identification of 16 microarray studies on metastatic melanoma published 2000 to 2011. Microarray data included in these studies are shown in Table S5. In the current study, we focused our attention on the differential gene expression between normal skin and/or benign nevi and metastatic melanoma. On this basis four microarray datasets were extracted (GEO access number: GSE7553, GSE4587, GSE4579, and GSE12391). An additional GSE22301 dataset was extracted from Rose et al. (2011), but while this study did not provide a gene signature of metastatic melanoma (and so was not included in the meta-analysis of 16 studies) it did include 14 samples of metastatic melanoma data and so was included in our integrative analysis. Thus, a total of five microarray datasets of normal and/or benign nevi and metastatic melanoma were used in this study (Table S6).

Genome-wide relative significance (GWRS) and Genome-wide global significance (GWGS) for integrated analysis of cross-laboratory microarray data

A relatively simple method of integrative meta-analysis was proposed by Rhodes et al. in 2002 that combines independent microarray studies based on the p-value of each individual gene: ${\displaystyle s_p=-\sum _{i=1}^{n}2log\left(p_i\right)}$ where p_i, i = 1–n, is the p-value of a gene in the i-th independent study. However, this method has at least two significant limitations: (1) many microarray studies are based on a small number of samples, for which the p-value can therefore be problematic and (2) the large variation in p-values across different studies leads to the data with smallest p-value determining the outcome of S_p.

We propose a new approach based on measuring the genome-wide relative significance (GWRS) and genome-wide global significance (GWGS) of expressed genes. We measure the GWRS of a gene using its ranking position (Jurman et al., 2008) on a genome-wide scale (r value) based on a differential expression measure, which can be the fold change, t-test p-value, SAM (Significance Analysis of Microarray data) p-value etc. Most existing meta analysis methods focus on the top-k genes (e.g. Jurman et al., 2008), while our method counts the ranking of genome-wide genes in total. Compared to the model of Rhodes and co-workers the proposed approach possess two important enhancements: (1) it can apply multiple different methods for measuring the degree of differential expression of a gene (e.g. fold change, t-test, Anova or SAM p-values) and (2) it uses a ranking r value instead of the test statistic (i.e., fold change, or p-value) to avoid the influence of high variation test statistics.

Data preparation

Pre-processing of microarray data is performed by extracting the expression value for each individual gene from the associated probe-sets. When a probe-set is mapped to multiple genes, e.g. ‘209994_s_at’ associated to two genes ‘ABCB1 / ABCB4’ in GSE4570, both genes are given the expression of the ‘209994_s_at’ probe-set.

For a gene appearing in multiple probe-sets, the most significant differential expressed probe-sets are assigned to this gene. We tested the results of using mean-, median-, and maxim-based methods to deal with the situations were multiple probe-sets are associated to a gene. We observed that the maxim-based method was able to retrieve the most significant probe-set of a gene, and would reflect our aim of extracting the most competitive genes across multiple studies. By contrast, use of a mean- or median-based probe-set value of a gene would drag the expression level down, and may introduce bias in follow-up analysis. As a result, a list of unique genes (G) from the datasets was retrieved. The number of datasets was denoted by n, while the number of unique genes across n datasets was denoted by m, i.e. m = |G|. The value ‘NA’ was applied in cases where a gene is absent from an individual study. We removed a gene from G where NA is bigger than δ (δ = 2 in this study), i.e. a gene was removed if it is absent for more than two of five datasets. This resulted in m = 24,097 and n = 5.

Measuring the GWRS of genes in each single microarray database

For each gene in the list of unique genes (G), we measured the degree of differential expression that can be measured by fold-change, t-test (p-value), SAM or other statistical test. However, fold-change is used in the current study, as our computational evaluation indicated that this produces more reliable results, probably due to the limited number of samples in some of the datasets. For each gene in G, we assigned a rank number (in descending order starting from 1 to m) according to their corresponding degree of differential expression i.e. a gene with a high degree of differential expression was ranked more highly and so with a smaller ranking number. An m∗n matrix (R) was thus created in which r_ij is the ranking number of the i-th gene in the j-th dataset. We measure the GWRS of the i-th gene in the j-th dataset by: ${\displaystyle s_{ij}=-2\hspace{0.17em}log\left(\frac{r_{ij}}{m}\right)}$ where r_ij, i = 1–m, j = 1–n, is the rank number of the i-th gene in the j-th study. The range of GWRS value (s_ij) is between 0 and −2log(1/m). For a gene with ‘NA’ value the s_ij is set to be ‘NA’.

Measuring the GWGS of a gene across multiple microarray datasets

We estimated the GWGS ($s_i^r$) of a gene based on its corresponding GWRS across n datasets, by ${\displaystyle s_i^r=\sum _{j=1}^n{\omega }_j{s}_{ij}}$ where ω_j represents the relative weight of the j-th dataset, and ${\sum }_{j=1}^n{\omega }_j=1$. The value of weight (ω_j) can be assigned based on the data quality of the j-th datasets (e.g. the level of data noise. The value of ω_j can also be used to reflect the differential importance of biopsy versus cell line samples that biological scientists may wish to take into account. In this study, we treated all the dataset equally, thus the weight of each datasets was set equally to be 1/n for j = 1–n. We also selected only the top 200 genes from the full gene list for further analysis (i.e. selected genes with the greatest s^r value) by empirical evaluation of the classification performance (accuracy ratio). This was determined using the ‘wrapper-feature selection’ after multiple rounds of gene addition (ranging from 20 genes up to 500 genes) in order to distinguish melanoma from normal skin/benign nevus. We observed that using more than 200 genes yielded no improvement in classification ratio values, and so we consider 200 genes as an optimal gene set with the smallest number of genes that still can achieve a similar level of classification performance.

Pathway analysis

We performed a pathway analysis to assess functional relevance of the new 200 gene signature based on the DAVID database (Hosack et al., 2003). DAVID provides a useful tool to analyze large gene lists, including via gene ontology and pathway analysis. We applied our top 200 genes to this database in order to detect potentially over-represented KEGG pathways. Before inputting into the DAVID database, we extracted the corresponding probe-sets of the 200 genes for the corresponding microarray platforms of each dataset. In comparison with the gene signature in the original 16 studies, we also extracted their associated probe-sets. We retrieved 31 pathways from the KEGG database where 12 genes (i.e., EGFR, FGFR2, FGFR3, IL8, PTPRF, TNC, CXCL13, COL11A1, CHP2, SHC4, PPP2R2C, and WNT4) in this 200-gene signature were found to closely interact with the 4 melanoma driver genes (see Results section).

Immunocytochemistry (ICC)

Primary epidermal melanocyte (EM) (female 44y), moderately pigmented human melanoma cells (FM55), and highly pigmented human melanoma cells (FM94) (melanoma cells were a gift of Dr Janis Ancans, University of Latvia) were cultured as previously described (Gledhill et al., 2010). The cells were fixed in ice-cold methanol (Sigma, Poole, Dorset, UK) for 10 min before air drying and rehydration in PBS. The cells were blocked with 10% donkey serum (DS) for 1 h, washed with PBS before incubation with respective primary antibodies to four test antigens from this 12-gene signature. These included: COL11A1 (Abcam, ab64883), CXCL13 (R & D Systems, AF801), PTPRF (NeuroMab, 75-193), SHC4 (Proteintech, 12641-1-AP), which were incubated overnight at 4 °C followed by secondary antibody (1:300) for 1 h (donkey anti-goat (Invitrogen, A11055), donkey anti-mouse (Invitrogen, A21202), donkey anti-rabbit (Invitrogen, A21206), Alexa green). The slides were cover-slipped by VECTASHIELD mounting medium with DAPI and photographed using a Nikon Eclipse 80i fluorescence microscope and imaged with a Nikon Digital Sight DS-U1 camera. A full assessment of all 12 proteins in our melanoma signature is beyond the scope of the current study, but will be assessed in detail in a follow-up studies.

Double immunohistochemistry (IHC)

Paraffin-embedded primary melanoma in situ (nose) and metastatic melanoma (lower leg) were deparaffinized and boiled in sodium citrate buffer (10 mM, 0.05% Tween 20, pH 6.0) for antigen retrieval. Acetone-fixed cryosections of normal human facial skin (Female 52 yrs) were used as control samples. All tissues were blocked with 10% donkey serum (DS) for 1 h, washed with PBS before 2 h incubation with NKi/beteb antibody raised against the melanocyte lineage-specific marker gp100 as a positive pigment cell control (Monosan; Mon7006-1) (1:15) followed by each of the 4 test antibodies at room temperature.

Data Access

The microarray data used in this study were retrieved from Gene Expression Omnibus (GEO) with the following access numbers: GSE4570, GSE4587, GSE7553, GSE12391, and GSE22301. The 16 signatures of melanoma reported in the literature between 2000 and 2011 were extracted from the associated publication and is presented in Table S2.

Gene signatures of melanoma (2000 to 2011) share few common genes

A meta-analysis conducted on gene signatures of metastatic melanoma reported in 16 independent microarray-based studies (ranging from 5 to 589 genes/study) from 2000 to 2011, showed remarkably few shared genes (Table 1, and Supplementary Information Table S1).

There were 84 genes common to two of the signatures (Scatolini et al., 2010; Jaeger et al., 2007), while 14 common genes appeared in three studies (Scatolini et al., 2010; Jaeger et al., 2007; Riker et al., 2008). Strikingly, while there were only 2 genes (KRT15, RORA) in common in four of the 16 studies (Scatolini et al., 2010; Jaeger et al., 2007; Riker et al., 2008; Smith, Hoek & Becker, 2005), we have recognized four genes in our 200 gene set (i.e. KRT15, MAGEA6, RORA and SULF1) that appeared in 4 different studies of the 16. No gene was common in five or more independent studies (Table S2). This finding suggested that there may be some fundamental issues with either the manner in which these microarray studies were designed, or with the meta-analyses conducted. On this basis we set about designing a new more robust model for meta-analysis.

Integrated analysis of cross-laboratory microarray data reveal a new melanoma gene signature

We applied our new approach to integratively analyze five independent microarray studies (Hoek et al., 2004; Smith, Hoek & Becker, 2005; Riker et al., 2008; Scatolini et al., 2010; Rose et al., 2011) (see Methods). The genome wide ‘global significance’ or GWGS of a gene (i.e., across all five datasets) was measured by the GWGS (s^r) as defined above (see Methods). A gene with a large s^r value is considered to be significant across multiple independent studies (i.e., globally significant). The 200 genes with largest s^r value were selected as the starting point for our new proposed gene signature of melanoma, as listed in Table 2 and Table S3. This set of 200 signature genes was empirically determined, based on the classification accuracy ratio after various rounds of gene additions (using the ‘wrapper feature selection’ approach) in order to distinguish melanoma from normal skin cells and/or benign nevus. As the classification accuracy ratio was improved very little by adding more than 200 genes, we applied this gene set as the smallest number of genes to retain the optimal classification accuracy performance.

Validation of a new 200-gene signature based on experimental studies reported in the literature

The 200 genes found to have genome-wide global significance in our study were compared with the gene signatures identified in previously-published reports (Table S5). Our new 200-gene signature was first validated by (i) comparing it with 16 signatures proposed in the referred to set of microarray studies (Table S1), (ii) checking if any experimental validation of these genes was published in the literature (PubMed, last access: 16 April 2012). This analysis revealed that (a) 85 genes in our 200-gene signature were reported in at least one of the 16 microarray studies, and (b) 21 genes of the 200-gene signature were reported in both microarray studies and wet-lab experimental studies (Table S4, labeled yellow background). We also found that 38 genes of this 200-gene signature were not reported in any of the 16 reference studies, but had in fact been previously validated in independent wet-lab studies (Table S4 and discussion section). Importantly, our new gene signature reported an additional subset of 77 genes that were not previously reported anywhere in the literature in association with melanoma (Fig. 1). The ranking positions of these 77 genes shows that 39% appear in the top 100 and 34% in bottom 50 (see Table S7). These genes may represent ‘novel genes’ as they were not previously identified in published microarray studies. We further investigated the characteristics of the 85 genes reported in at least 1 of the 16 reference microarray studies (Table S3). Forty-four were reported in ≥ 2 studies, while 17 genes have been reported in ≥ 3 of the 16 studies (Table S3). KRT15, MAGEA6, RORA and SULF1 were the most frequently reported genes appearing in 4 of the 16 studies. Thus, using our method, we are able to pick up 4 of the 7 most frequently reported genes in the 16 studies by using just our top 200 genes (i.e., 30% less than the next best list of 308 genes in Jaeger et al. (2007)). In this way the methodology to select the top 200 genes in our study is more powerful than previously reported on the component 16 published signatures used for the source data (Table 1).

Interaction of a new 200-gene signature with melanoma ‘driver’ genes informs a new signaling network in melanoma

We investigated the interaction between genes within our 200-gene signature with the four known melanoma ‘driver’ genes (i.e., NRAS, BRAF, MITF and cKIT). Of these driver genes, NRAS is mutated in 13–25% of melanoma cases (Goel et al., 2006; Schubbert, Shannon & Bollag, 2007), while BRAF (located downstream of NRAS), is mutated in up to 45% of malignant melanomas (Hocker & Tsao, 2007; Flaherty & McArthur, 2010). MITF, a master transcription factor in melanocyte function, cooperates when mutated with BRAF in melanomagenesis (Garraway et al., 2005; Taylor et al., 2011). Recent studies show that mutant cKIT can activate the Ras/Raf/Mek/Erk pathway and also activate MITF (Monsel et al., 2009; Phung et al., 2011). The four well-known melanoma driver genes did not appear on our list. This is due most likely to these four driver genes being associated with melanoma at the gene mutation level, rather than at the gene expression level.

We retrieved 31 pathways from the KEGG database where 12 genes in our proposed 200-gene signature were found to closely interact with the 4 melanoma driver genes in the MAPK, Ca²⁺ and WNT signaling pathways (Table 3). These 12 genes are EGFR, FGFR2, FGFR3, IL8, PTPRF, TNC, CXCL13, COL11A1, CHP2, SHC4, PPP2R2C, and WNT4. Based on these interactions we propose a new signaling network for melanoma (Fig. 2). Of these 12 genes, CXCL13, SHC4, WNT4 and CHP2 were detected only using our computational method (i.e., not reported before in melanoma microarray studies) but exhibit important positions in melanoma driver gene signaling pathways (Fig. 2). The biological pathways involving chemokine receptors, WNT, Ca²⁺ and MAPK signaling will have implications for melanomagenesis and metastatic progression.

Experimental validation of a MAPK pathway-associated subset in our 12-gene melanoma signature

Four genes in our proposed 12-gene biomarker signature that appear in the MAPK signaling pathway (i.e., COL11A1, CXCL13, PTPRF, and SHC4) were selected for laboratory validation. Note that COL11A1, CXCL13, and PTPRF have not previously been reported to be associated with melanoma experimentally. COL11A1, CXCL13, PTPRF, and SHC4 were found to be over-expressed in two human melanoma cell lines (i.e., FM55 and FM94) compared to normal human epidermal melanocytes in vitro (Fig. 3). A significant degree of heterogeneity was observed in the expression pattern for these markers. For example, COL11A1, a secreted collagen protein, was observed at low levels in the cytoplasm of normal melanocytes, but more intensely in the perikayon of moderately-pigmented FM55 melanoma cells, and unexpectedly exhibited a nuclear/nuclear membrane association in the pigmented FM94 melanoma cells. Similarly, a weak cytoplasmic localization of CXCL13 in normal melanocytes appeared to shift towards the perikayon and nucleus of FM55 and FM94 melanoma cells respectively, as evidenced by co-localization with DAPI staining. Low level PTPRF expression in normal epidermal contrasted with higher expression (both cytoplasmic and nuclear) in melanoma cells. Finally, SHC4 expression was membranous in normal melanocytes contrasting with some punctuate nuclear membrane expression in melanoma cells (Fig. 3).

The expression of these four proteins was also assessed in normal human healthy skin and in melanoma patient tissue (both primary and metastatic melanoma). Using double immunofluorescence with a melanocyte lineage marker gp100, we assessed the relationship of the four test proteins with melanocytes or melanoma cells in these tumor biopsies. We included primary melanoma (in addition to metastatic melanoma) in our immunohistochemistry validation study because the expression levels for the 12 genes in our signature exhibited several fold level changes between primary melanoma and normal skin/benign nevi across 5 microarray datasets (Table S8).

COL11A1, CXCL13 and PTPRF were not detected in normal human epidermal melanocytes in situ (Fig. 4a). Some low level expression of SHC4 was detected in these normal pigment cells. By contrast, COL11A1 was expressed intensely by melanoma cells located in the dermis of both primary and metastatic melanoma (Fig. 4b). CXCL13 was strongly expressed in a minor subpopulation of tumor cells in primary melanoma, while a greater fraction of cells in metastatic melanoma tissue expressed this protein. By contrast, PTPRF was intensely expressed in the majority of tumor cells of both primary and metastatic melanoma cells. Finally SHC4 was found to be expressed in minor fraction of primary gp100-positive melanoma, but in most metastatic gp100-positive melanoma cells.

Computational evaluation of the robustness of a proposed 12-gene biomarker signature in distinguishing melanoma from normal skin and/or benign nevi

A computational evaluation of robustness of the proposed 12-gene signature, based on melanoma driver gene association, was performed for distinguishing melanoma from normal skin and/or benign nevi using cross-laboratory published data. This data evaluation is important to verify the robustness of a new biomarker for potential diagnostic application and/or possible therapeutic development. The support vector machine (so-called SVM model) classification model (Brown et al., 2000) and the ‘leave-one-out method’ are used to classify microarray datasets (Hoek et al., 2004; Smith, Hoek & Becker, 2005; Riker et al., 2008; Scatolini et al., 2010; Rose et al., 2011). Our results showed that these 12 genes achieved excellent classification accuracy ratios across these five datasets (i.e., average of 99.1%, Table 4). This result indicated that our 12-gene biomarker achieved a classification accuracy ratios that was identical or near identical to the classification accuracy ratios of the original individual studies. Importantly, the 12-gene biomarker signature achieved a much better performance on average than the signatures of Smith, Hoek & Becker (2005), Riker et al. (2008) and Scatolini et al. (2010), and very slightly less (0.44% less) classification accuracy than the signature of Hoek et al. (2004). It should be noted that the signature of Hoek et al. (2004) consisted of 589 genes, while our biomarker signature is very much shorter at just 12 genes.