Epigenomic and genomic analysis of transcriptome modulation in skin cutaneous melanoma

DNA MET and copy number dysregulation at the transcriptomic level

MET, CNV, and mRNA expression profiles were obtained from 466 SKCM samples in TCGA. Original data were then preprocessed as described in Materials and Methods. Correlation coefficients between DNA MET or CNV profiles and mRNA expression profiles were calculated to assess the effects of epigenomic and/or genomic aberrations. Correlation coefficient r values were normalized using Fisher’s Z-transformation to stabilize variance.

Consistent with our previous findings, overall correlation coefficients between DNA CNVs and expression profiles showed a marked right-sided skew (skewness = 0.67425, p < 1e^-5). By contrast, correlation coefficients between DNA MET patterns and expression profiles displayed a left-sided skew (skewness = -0.4274, p < 1e^-5) (Figure 1A), indicating that abnormalities in DNA MET and CNVs negatively and positively regulated transcription, respectively.

Figure 1. Identification of DNA copy-number-correlated (CNVcor) and DNA methylation-correlated (METcor) genes in SKCM. (A) Distribution of correlation coefficients between mRNA expression levels and DNA copy number or DNA methylation across samples. (B) Venn diagram showing counts of CNVcor genes, METcor genes, and overlapping genes. (C) Proportional frequencies of CNVcor genes against total gene counts in each chromosome arm. (D) Proportional frequencies of METcor genes against total gene counts in each chromosome arm. (E, F) Genomic positions of DNA methylation probes are categorized based on positional relations with CpG islands (right) and genes (left).

The METcor and CNVcor gene sets included many individual genes (Supplementary Tables 2–3); only genes with notable correlations with overall survival (OS) were included in subsequent analysis (log rank p < 0.05). Positively-correlated gene signatures were typically included in the DNA copy number gene set (CNVcor, n = 385), while negatively-correlated gene signatures were typically included in the DNA MET gene set (METcor, n = 892). CNVcor genes were associated with CNV-dependent dysregulation at the transcriptional level, whereas METcor genes were associated with MET-dependent transcriptional dysregulation at the epigenetic level. In addition to genes that belonged to either the METcor or the CNVcor set, 216 overlapping genes were included in both sets, implying that both METcor and CNVcor genes were specifically associated with dysregulation at the transcription level (Figure 1B).

As was found in a previous study, chromosome 1 was rich in CNVcor genes (Figure 1C), suggestive of regional sensitivity of gene expression to DNA dosage [20, 21]. METcor genes were also preferentially located in certain chromosomal regions, such as chromosomes 17 and 1 (Figure 1D and Supplementary Figure 1); in addition, most METcor genes were protein coding genes (MET on the promoter region) (Figure 1E) and located primarily within CpG islands (Figure 1F). These results indicate that METcor and CNVcor genes were important contributors to transcriptional dysregulation in SKCM.

Molecular subtypes based on METcor and CNVcor genes

Next, we examined whether METcor and CNVcor gene expression could predict prognostic subgroups. Each gene set profile was examined using NMF clustering analysis, with the cluster number k set at 2-10; k values were then determined for all profiles (k=3 for MET and k = 2 for CNV) (Figure 2A, 2B and Supplementary Figures 2–5). Surprisingly, subtype identifications based on CNVcor genes overlapped to a large extent with those based on METcor genes (p < 1e^-5, χ²-test), suggesting that METcor and CNVcor genes both contribute to regulation of SKCM (Figure 2E, 2F). Moreover, Kaplan-Meier (KM) curve analysis indicated that subtypes identified based on either METcor or CNVcor genes were associated with patient OS (Figure 2C, 2D, p < 0.05).

Figure 2. Identification of SKCM molecular subtypes using CNVcor and METcor genes. (A, B) Plots show the non-negative factorization (NMF) cluster results for CNVcor genes in CNV data (A) and for METcor genes in MET data (B). (C, D) Kaplan–Meier plot analyses of differences in OS among subtypes identified by NMF clustering of CNVcor (C) and METcor (D) genes. (E, F) Subtypes based on CNVcor genes overlapped to a great extent with those based on METcor genes.

The comprehensive clustering method iCluster was adopted to integrate genomic information regarding DNA MET and CNVs as well as mRNA expression. Clustering analysis was performed using a cluster number k of 2/3. Twenty iterations of clustering analysis at K=2 (category 3) and at K=3 (category 4) were performed to generate optimal iCluster clustering results. Our results indicated that K=3 generated more stable clustering results than K=2 (Supplementary Figures 6, 7). Samples were therefore clustered into four subclasses: iC1-iC4 (n=99, 108, 113, and 146, respectively). Clustering results for these four subclasses are displayed in Figure 3A, 3B; clustering results for each sample are presented in Supplementary Table 4.

Figure 3. Expression patterns of subtypes based on CNVcor and METcor genes. Subtypes based on CNVcor (A) or METcor (B) genes using NMF cluster methods are indicated with colored bars.

According to KM curve analysis, OS was poorest in the iC3 subtype (p < 0.05, Figure 4A). Comparisons of patient OS among the four subgroups (Figure 4B and Supplementary Figure 8) indicated that the iC1 and iC3 subgroups differed most in terms of prognosis (p < 0.001). In addition, the subtypes identified using iCluster largely overlapped with those identified based on METcor and CNVcor genes (p < 1e-5, χ²-test, Figure 4C, 4D). These results suggest that integrated analysis of METcor and CNVcor genes identified clinically-relevant molecular subtypes in which epigenomic and genomic transcriptional dysregulation influences prognosis.

Figure 4. Identification of SKCM molecular subtypes using iCluster analysis. (A) Kaplan–Meier plot analyses of differences in OS among subtypes identified by iCluster (iC1, iC2, iC3, and iC4). (B) Kaplan–Meier plot analyses of OS in iC3 and iC1 subtypes. (C, D) The subtypes determined through iCluster analysis overlapped extensively with those based on CNVcor (C) or METcor (D) genes.

Abnormalities in DNA MET and CNVs are linked

After correcting for batch effects, occurrence rates of genome-wide abnormalities in DNA MET and CNVs were compared. DNA copy-number loss (CNVloss, β<-0.3) and gain (CNVgain, β>0.3), as well as DNA hypomethylation (METhypo, β<0.2) and hypermethylation (METhyper, β>0.8), were also identified based on the predetermined 0.3-fold change threshold and were compared with the average values of each probe in normal tissue. Our results (Supplementary Table 5) indicate that CNVgain and CNVloss frequencies were strongly correlated (p < 1e-5, Figure 5A). In addition, METhyper and METhypo frequencies were also correlated (p < 1e-5, Figure 5F). Abnormalities in directional CNVloss, CNVgain, METhypo, and METhyper were tightly correlated, indicating that each correlation was independent of directional abnormality (Figure 5B–5E, p < 0.001). In summary, our findings indicated that SKCM cases with higher DNA CNV frequencies also tended to have higher frequencies of abnormal DNA MET.

Figure 5. Abnormal DNA copy numbers and DNA methylation are linked in SKCM. (A, F) Abnormalities in DNA copy numbers or DNA methylation were determined based on a cutoff fold difference >0.2 compared to average values in non-tumor tissues. DNA copy-number gains (CNVgain) and losses (CNVloss) and DNA hypermethylation (METhyper) and hypomethylation (METhypo) are shown for each sample. (B–E) Plots show the pairwise frequencies of CNVgain, CNVloss, METhyper, and METhypo genes in individual samples.

Identification of key molecular features of SKCM subtypes

Differences in clinical characteristics (such as Stage, TNM, Primary Site, and Gender) were examined among the four subtypes identified above. There were no statistically significant differences in clinical characteristics among the four SKCM subtypes (Supplementary Figure 9). In addition, the tumor immune estimation resource (TIMER) approach was used to analyze differences in samples belonging to the four subtypes (Supplementary Table 6) [22]. Six immune cell scores for samples in the iC1 subtype that had the best prognosis were markedly elevated compared to scores for samples in other subtypes (Figure 6A, 6B, p < 0.01). These results suggest that degree of immune cell infiltration or immune microenvironment in SKCM are correlated with DNA MET or CNV levels.

Figure 6. Identification of key immune features for SKCM subtypes. Immune scores were calculated for samples from the four subtypes and were compared using the TIMER (tumor immune estimation resource) approach. Scores were determined for six immunocytes in each sample (A); these scores were compared among subtypes (B).

Next, we examined heterogeneity in DNA MET, CNVs, and gene expression among iC3 and iC1 subtype samples. Three categories each were created for DNA CNV and MET levels: Normal, Loss, and Gain for CNV and Normal, HypoMethy, and HyperMethy for MET. Fisher’s exact test was used to identify DNA MET or CNV genes with levels that differed between the iC3 and iC1 subtypes [23]; the results are shown in Supplementary Tables 7, 8. Differentially expressed genes (DEGs) between the iC3 and iC1 subtypes were identified using DESeq2 (Supplementary Table 9, p < 0.05) [24]. MET, CNV, and expression results for genes with marked differences among samples from the four subtypes are shown in Figure 7. A total of 221 genes with different MET, CNV, and expression levels between the iC3 and iC1 subtypes were examined in univariate survival analysis to identify differences in prognostic characteristics among the subtypes. The results indicated that 146 of these genes were correlated with prognosis (Supplementary Table 10, log rank p < 0.05). Among them, increased hypermethylation and decreased expression of the GBP5, CD8A, SAMHD1, and KIAA0040 genes were associated with poorer outcomes in the iC3 subtype compared to the iC1 subtype. Samples were then assigned to low (L1), medium (L2), or high (L3) groups based on the expression of KIAA0040, SAMHD1, CD8A, and GBP5. The results indicated that prognostic outcomes were positively correlated with the expression of CD8A, GBP5, KIAA0040, and SAMHD1 (Figure 8A–8D), suggesting that CD8A, GBP5, KIAA0040, and SAMHD1 expression levels were associated with DNA CNV or MET level. Next, associations between the abovementioned 146 genes and patient prognosis were evaluated in the GEO GSE65904 SKCM dataset (Supplementary Table 10). Seventy-eight genes were correlated with prognosis in that dataset; survival curves for the top 20 genes are shown in Supplementary Table 11 and Supplementary Figure 10.

Figure 7. Identification of general molecular features of SKCM subtypes. (A) Distribution of DNA CNVs in iCluster subtypes. (B) Distribution of DNA MET in iCluster subtypes. (C) Heatmap of DEGs among iCluster subtypes.

Figure 8. Identification of key molecular features in SKCM subtypes. (A) GBP5 MET, CNV, and expression levels in the iC1 and iC3 subtypes are shown in the left panel. Samples were divided among high (L3), medium (L2), and low (L1) groups based on GBP5 gene expression; Kaplan–Meier plot analysis for these groups is shown in the right panel. (B) CD8A MET, CNV, and expression levels in the iC1 and iC3 subtypes are shown in the left panel; Kaplan–Meier plot analysis for L1, L2, and L3 CD8A groups is shown in the right panel. (C) SAMHD1 MET, CNV, and expression levels in the iC3 and iC4 subtypes are shown in the left panel; Kaplan–Meier plot analysis for L1, L2, and L3 SAMHD1 groups is shown in the right panel. (D) KIAA0040 MET, CNV, and expression levels in the iC3 and iC4 subtypes are shown in the left panel; Kaplan–Meier plot analysis for L1, L2, and L3 KIAA0040 groups is shown in the right panel.

Finally, SKCM mutation profiling data were examined to explore associations with the subclassifications identified in this study. Synonymous mutations were removed, and both nonsense and missense gene mutations were included. Overall, mutation frequencies differed significantly among the various subtypes. A selection of 85 genes with mutation frequencies that differed markedly (P < 0.01) between the iC1 and iC3 subtypes based on Fisher’s test is shown in Figure 9; mutation frequency data are shown in Supplementary Table 12, and Fisher’s test results are shown in Supplementary Table 13. For example FAM135B gene mutation frequency was much higher in the iC3 subtype than in the iC1 subtype (p < 0.01). Interestingly, FAM135B gene mutation frequency was also increased in the iC2 and iC4 subtypes, which had poor prognoses, compared to the iC1 subtype (p < 0.05). In summary, these results suggest that DNA copy number- and MET-related molecular subtypes of SKCM are associated with differences in FAM135B gene mutation frequency, which might regulate progression of SKCM.

Figure 9. Differentially mutated genes among the SKCM subtypes.

DNA copy number, DNA MET, and mRNA expression profiles

The SKCM dataset was downloaded from TCGA. DNA copy number, DNA MET, mRNA expression patterns, and mutation data for SKCM samples were extracted through the official TCGA data portal. A total of 466 samples for which matched DNA copy number, DNA MET, and mRNA expression pattern data were available were included in this study (Supplementary Table 1). Gene expression profiles were normalized using quantile normalization and log2 transformation, and were later aggregated according to official HUGO symbols. Each expression profile was then normalized to mean expression of the probe gene in normal tissue so that it represented the fold change in cancer tissue relative to normal tissue. The circular binary segmentation algorithm from the R package “DNAcopy” library was then used to map genetic DNA copy numbers to CNV data in every sample according to the segmented CNV data [32]. For DNA MET patterns, probe β-values were filtered to remove probes located on sex chromosomes. Then, probes located in CpG island-associated regions were mapped to related genes, including Shelf, CpG islands, Shore regions, first-exon regions, differentially methylated regions, 5’-UTRs, and gene promoter regions containing 2500 upstream bases from TSS. Probes with >30% missing values across samples were removed from each processed profile.

The R package “liftOver” library was used to align probe genomic coordinates to the hg38 human reference genome in each dataset. Probes were then matched to corresponding probes in the mRNA expression profile data. In addition, cancer-specific alterations were identified by subtracting the average intensity of the probe in normal tissue. After probes with >50% missing values and those located on sex chromosomes had been removed, data were analyzed using the sk-nearest neighbor algorithm. Pairwise Pearson’s correlation coefficients were then calculated for each gene within the paired profiling patterns of CNV and expression and MET and EXP. If multiple probes had been mapped to the same gene, the probe with the mean or lowest correlation coefficient was used as the representative pair-matched probe for MET and CNV profiles, respectively.

Clustering analysis of genome patterns from different aspects

Stable sample clusters were identified through negative matrix factorization (NMF) clustering analysis according to the “brunet” method and 50 iterations using METcor and CNVcor genes. The cluster number k was set between 2 and 10, and the optimal cluster number was computed based on the monitored consensus map as well as cluster cophenetic correlations. In addition, average silhouette width was determined for the consensus membership matrix using the R package “NMF.” For each member, the smallest cluster number was set to be 10. The “iCluster” R package was used to perform comprehensive clustering analyses of DNA CNV, MET, and mRNA expression profiles using default parameters and 20 iterations.