Exome sequencing of desmoplastic melanoma identifies recurrent NFKBIE promoter mutations and diverse activating mutations in the MAPK pathway.

Desmoplastic melanoma is an uncommon variant of melanoma with sarcomatous histology, distinct clinical behavior and unknown pathogenesis. We performed low-coverage genome and high-coverage exome sequencing of 20 desmoplastic melanomas, followed by targeted sequencing of 293 genes in a validation cohort of 42 cases. A high mutation burden (median of 62 mutations/Mb) ranked desmoplastic melanoma among the most highly mutated cancers. Mutation patterns strongly implicate ultraviolet radiation as the dominant mutagen, indicating a superficially located cell of origin. Newly identified alterations included recurrent promoter mutations of NFKBIE, encoding NF-κB inhibitor ɛ (IκBɛ), in 14.5% of samples. Common oncogenic mutations in melanomas, in particular in BRAF (encoding p.Val600Glu) and NRAS (encoding p.Gln61Lys or p.Gln61Arg), were absent. Instead, other genetic alterations known to activate the MAPK and PI3K signaling cascades were identified in 73% of samples, affecting NF1, CBL, ERBB2, MAP2K1, MAP3K1, BRAF, EGFR, PTPN11, MET, RAC1, SOS2, NRAS and PIK3CA, some of which are candidates for targeted therapies.

Desmoplastic melanoma is an uncommon variant of melanoma with sarcomatous histology, distinct clinical behavior, and unknown pathogenesis[1-3]. We performed low-coverage genome and high-coverage exome sequencing of 20 desmoplastic melanomas, followed by targeted sequencing of 293 genes on a validation cohort of 42 cases. A high mutation burden (median 62 mutations/Mb) ranked desmoplastic melanoma among the most highly mutated cancers[4]. Mutation patterns strongly implicate UV-radiation as the dominant mutagen[5], indicating a superficially located cell of origin. Novel alterations included recurrent promoter mutations of NF-kappa B inhibitor epsilon, NFKBIE (IkBε), in 14.5% of samples. Commonly mutated oncogenes in melanomas, in particular BRAF and NRAS, were absent. Instead, other genetic alterations known to activate the MAPK and PI3K signaling cascades were identified in 73% of samples, affecting NF1, CBL, ERBB2, MAP2K1, MAP3K1, BRAF, EGFR, PTPN11, MET, RAC1, SOS2, NRAS, and PIK3CA, some of which are candidates for targeted therapies.

Desmoplastic melanomas comprise 4% of all primary melanomas[1]. They most often occur in chronically sun-exposed skin of older individuals. Clinically, they usually present as unpigmented scar-like indurations, delaying their detection. Histologically they are primarily dermal-based tumors composed of a paucicellular proliferation of spindle-shaped cells situated within an abundant desmoplastic stroma and can be easily misdiagnosed. Some desmoplastic melanomas grow extensively along or within nerves, a feature termed neurotropism. Primary melanomas composed predominantly (>90%) of desmoplastic melanoma are classified as “pure desmoplastic melanoma", whereas primary melanomas that exhibit areas characteristic of desmoplastic melanoma amounting to <90% of the tumors (the remainder being composed of non-desmoplastic melanoma), are classified as “mixed desmoplastic melanoma”. Patients with the pure subtype of desmoplastic melanoma have a lower rate of nodal metastasis and better survival[2,3,6,7].

Despite the fact that desmoplastic melanoma is a deadly and commonly misdiagnosed cancer, there are no known genetic drivers[2]. Several studies have failed to identify common pathogenic mutations[8,9]. This may be due to the small number of cases analyzed to date as well as the technical limitations of studying tumors like desmoplastic melanoma that exhibit a low ratio of tumor to stromal cells. We sought to characterize the genome-wide mutational landscape of desmoplastic melanoma in order to identify genetic alterations that underlie their unique biology and serve as potential diagnostic biomarkers or therapeutic targets.

Low-coverage genome sequencing (13×) and high-coverage exome sequencing (89×) were performed on a discovery set of 20 fresh-frozen tumors and matched normal DNA (Supplementary Table 1). In addition, we sequenced neoplastic and non-neoplastic tissue from 42 formalin-fixed paraffin-embedded primary desmoplastic melanomas as a validation cohort. Targeted sequencing of 293 genes (216×, Supplementary Table 2 for gene list) that included the top candidates nominated from the discovery set was performed on the validation cohort. High-resolution aCGH (180K-1M feature arrays) was performed on samples from both cohorts. Using these approaches, we were able to determine point mutation and copy number information in 62 desmoplastic melanomas.

The median number of mutations per megabase was 62 (Fig 1), ranking desmoplastic melanoma among cancers with the highest known mutation burdens. This contrasts with “conventional” cutaneous melanomas (approximately 15 mutations per megabase[10,11]) and the majority of solid cancers (approximately 2 mutations per megabase[4,5]). In desmoplastic melanoma 88% of mutations were C>T transitions and favored di-pyrimidines, implicating UV-radiation as the dominant mutagen (Fig 1). There were some notable exceptions to these patterns. Two tumors arose from sun-shielded sites and had the lowest mutation burdens; one of these occurred in a patient in whom we identified a germline CDKN2A mutation. Tumors arising in younger patients tended to have lower mutation burdens (Fig 1) (p =2×10-3, t-test). Pure and mixed desmoplastic melanomas were genetically similar.

Figure 1

Desmoplastic melanomas have a substantial point mutation burden consistent with UV-radiation induced damage

The 62 tumors are ordered by their mutation burden (top panel) with the mutation types annotated (bottom panel). Dashed line corresponds to the mutation burden observed in sun-exposed non-desmoplastic melanoma: 15 mutations/Mb[10,11]. In the bottom panel, C>T mutations following a purine (*) or a pyrimidine (**) are distinguished. Tumors from patients older than 55 years of age had significantly more mutations as compared to those from younger patients (p=2×10-3, t-test). Cases from the discovery cohort are marked with a “D”.

Overall, desmoplastic melanomas had fewer copy number alterations (CNAs) than other melanoma subtypes (Supplementary Figure 1). In spite of the low overall copy number burden, several samples had focal copy number alterations (Supplementary Figures 2-3). Focal amplifications affected the following genes: EGFR, CDK4, MDM2, TERT, MAP3K1 (3 cases each); MET, YAP1, NFKBIE (2 cases each); CCND1, MYC, and SOS2 (one case each). Immunohistochemistry for EGFR, CDK4, MDM2, MET, YAP1, and CCND1 confirmed increased protein expression in select cases (Supplementary Figure 4). Focal deletions affected CDKN2A (11 cases) and NF1 (4 cases). Loss of p16 expression was confirmed by immunohistochemistry in select cases (Supplementary Figure 5).

The high mutation burden made it difficult to identify potential driver mutations among the numerous somatic mutations. We searched for recurrent mutations clustering at specific base pairs (hotspots) and genes with a disproportionately high frequency of loss-of-function mutations to identify true driver mutations[10-12].

Several tumors harbored the following oncogenic hotspot mutations known to occur in other cancers: ERBB2 (n=4); MAP2K1, PTPN11, PPP6C, and RAC1 (n=2); EZH2, IDH1, PIK3CA, NRAS, BRAF, BRAF, and BRAF (n=1) (Supplementary Tables 3-4). TERT promoter mutations were only recently discovered[13,14]; the relevant region could only be sequenced in a subset of samples, in which 85% (17/20) harbored a mutation (Supplementary Figure 6). To identify novel pathogenic mutations, we searched for highly recurrent mutational hotspots that have not been previously characterized. The most recurrent mutational hotspot, observed 9 times, affected the gene NFKBIE (Fig 2A). NFKBIE also harbored recurrent mutations at another nearby position in two additional instances (Fig 2A).

Figure 2

Nomination of driver mutations in desmoplastic melanoma

A. The number of specific mutations detected (y-axis) is stratified by those mutations' occurrences across samples (x-axis). For example, 3,364 mutations occurred only once, whereas there were 151 mutations occurring in two samples and a single mutation was observed 9 times, affecting the NFKBIE genetic locus. There was also a secondary hotspot in NFKBIE, 15 basepairs from the more common mutation site, which was mutated in two samples. B. Q-Q plot of loss-of-function burdens compared to expected loss-of-function burdens, calculated as described. Solid and dotted lines correspond to false discovery rates (FDRs) of 1.0 and 0.5, respectively. The most significant genes are labeled. C. Tumor suppressor candidates have an increased proportion of damaging mutations and fully clonal mutant allele frequencies (MAFs), undergoing loss of heterozygosity in some cases. Left panel: Fraction of mutation categories compared to all mutations. Right panel: Normalized MAFs (calculated as described) of candidate mutations compared to all mutations. Red vertical bars indicate average MAFs.

To identify tumor suppressor candidates, we looked for genes that were enriched for loss-of-function mutations, as described in the Methods. Briefly, truncating mutations (nonsense, splice-site, or frame-shift) and missense mutations predicted to be damaging were nominated as candidate driver mutations. Mutations that had undergone loss of heterozygosity were particularly scrutinized. Taking these criteria into account, we determined genes for which the burden of loss-of-function (LOF) exceeded what would be expected by chance by comparing it to LOF burdens generated from permuted data (Fig 2B). The genes implicated by this approach were TP53, CDKN2A, NF1, ARID2, FSIP1, CBL, FBXW7, CDH2, IL36A, PAK3, RB1, and ARID1A (Fig 2C). Using immunohistochemistry, we confirmed that representative mutations affecting TP53, CDKN2A, RB1, and ARID1A resulted in concurrent alterations at the protein level (Supplementary Figure 5).

Our analysis implicated several cancer genes unique or enriched in desmoplastic melanoma. CBL is an E3 ubiquitin ligase that targets several receptor tyrosine kinases (RTKs) for degradation, and its loss is associated with increased RTK signaling[15]. Germline CBL mutations affecting its zinc finger domain are associated with a variety of cancer and developmental disorders, such as Noonan syndrome[16]. Noonan syndrome can be caused by mutations in other genes, including PTPN11, NRAS, KRAS, and BRAF, all of which have been associated with melanocytic neoplasia[17]. In desmoplastic melanoma, CBL harbored frequent truncating and damaging missense mutations in the absence of any synonymous mutations (Fig 3A), indicating that it may act as a tumor suppressor gene. CBL mutations are infrequent in published melanoma exome sequencing studies. This could be because non-desmoplastic melanomas generally harbor BRAF and NRAS[10,11] mutations, which functionally reside in the same signaling pathway as CBL.

Figure 3

Genetic alterations of CBL, MAP3K1, FBXW7, and NFKBIE

A. Enrichment of damaging mutations of CBL in desmoplastic melanoma. B. Recurrent amplification of MAP3K1. C. Enrichment of damaging mutations of FBXW7 in desmoplastic melanoma. D. Recurrent amplification of NFKBIE.

Another gene implicated in the MAPK pathway was MAP3K1[18,19]. We observed highly focal amplifications of MAP3K1 in three tumors (Fig 3B). Although genetic alterations of MAP3K1 have not been reported in melanoma, a Sleeping Beauty screen in a melanoma mouse model driven by BRAF identified MAP3K1 as a melanoma oncogene[20]. In that study, tumors with MAP3K1 insertions arose exclusively in melanocytes that failed to activate the conditional BRAF allele, or they arose in BRAF control mice, indicating that MAP3K1 activation can substitute for BRAF activation. MAP3K1 amplifications may therefore represent an equivalent driver mutation in desmoplastic melanoma.

FBXW7 is an E3-ubiquitin ligase responsible for MYC and CCNE1 degradation[21,22]. Truncating or damaging missense mutations striking the critical WD domains of FBXW7 are common in several cancers[23-25]. In our cohort, 11% of tumors harbored nonsense or damaging missense mutations, often involving the WD domains, while synonymous and conservative missense mutations were absent (Fig 3C). FBXW7 mutations are also present, albeit at a much lower frequency, in published melanoma exome sequencing studies [10,26-29]. Therefore, FBXW7 is a tumor suppressor across multiple melanoma subtypes with an increased frequency in desmoplastic melanoma.

NFKBIE (IkBε) inhibits downstream Nuclear Factor kappa B (NFκB) signaling by sequestering NFκB transcription factors in the cytoplasm[30]. The NFKBIE locus was focally amplified in two samples (Fig 3D), and showed recurrent mutations at several nearby hotspots, implicating NFKBIE as a candidate oncogene. The clustered mutations are not found in COSMIC nor TCGA studies of any cancer, including sequencing studies of melanomas mostly from intermittently sun-damaged skin. We sequenced NFKBIE in an extension cohort of diverse melanomas and found 6 non-desmoplastic melanomas with NFKBIE mutations (see supplementary text). Similar to desmoplastic melanomas, the NFKBIE-mutant non-desmoplastic melanomas did not harbor BRAF or NRAS mutations, and had a very high mutation burden with evidence of UV-radiation induced mutational damage.

Overall, we found 20 clustered mutations in 15 tumors (Fig 4). Mutations from all cases were validated using at least two sequencing assays (Supplementary Table 6 and Supplementary Figures S7-S9). Remarkably, five of the 15 tumors had two mutations each, and in all five cases the mutations affected opposing alleles (Fig 4 and S8). Four of the mutations from the 10 tumors with only a single NFKBIE mutation had undergone loss of heterozygosity, as evidenced by elevated frequencies of the mutant allele (Fig 4, Supplementary Table 6). Collectively, these results strongly suggest selection for bi-allelic NFKBIE mutations.

Figure 4

Recurrent mutations affect the promoter of NFKBIE

Mutations are annotated over the entire NFKBIE genetic locus (top panel) and a zoomed inset (bottom panel) with selected tracks featured as described in the main text.

There are two gene models for NFKBIE. The main mutational hotspot resided over the coding region of the long isoform and the promoter of the short isoform (Fig 4, NFKBIE genes track). We mined Illumina Human Bodymap 2.0 RNA-Seq data to determine the tissue distribution of NFKBIE isoforms. The short isoform of NFKBIE is ubiquitously expressed in all tissues, while the long form is restricted to brain tissue (Fig 4, NFKBIE genes track). Similarly, RNASeq data from the melanoma TCGA project exclusively showed expression of the short isoform (Fig 4, Melanoma TCGA RNA-Seq track). We identified two non-desmoplastic melanoma cell lines that harbored NFKBIE hotspot mutations (M257 and M375, Fig 4 Mutations track). Using isoform-specific RT-PCR (Supplementary Figure 10A) and Western blot (Supplementary Figure 10B-C) we found that these cell lines also only expressed the short isoform. In conjunction, the mutational hotspot of NFKBIE affects the promoter region of the relevant isoform (Fig 4, Promoter track) in an area that is highly conserved across multiple species (Fig 4, Conservation tracks) and predicted to affect binding sites for 32 transcription factors, including the consensus binding motifs for GABPA and ELF1 (Fig 4, Transcription Factor tracks).

Many melanoma cell lines have been reported to show nuclear localization of NFκB transcription factors, suggesting that NFκB signaling is active in melanoma[31-34]. By contrast we found that NFκB nuclear translocation was absent in the two NFKBIE-mutant cell lines, M375 and M257 (Supplementary Figure 10D). This result is consistent with the proposed gain-of-function role for NFKBIE resulting in inactivation of NFκB signaling. However, future studies will be necessary to dissect the specific mechanism by which these promoter mutations modulate signaling.

Many of the individually nominated driver genes are components of critical signaling pathways in melanoma and cancer in general. The RTK→Ras→MAPK, and →PI3K signaling cascades harbored genetic alterations predicted to lead to its activation in 73% of tumors (Fig 5). Concordantly, immunohistochemistry for phospho-ERK was ubiquitously positive in desmoplastic melanomas (Supplementary Figure 11), confirming the importance of MAP-kinase pathway activation in this melanoma type as well. In spite of this similarity with other melanoma types, the genes affected within the MAP-kinase pathway in desmoplastic melanoma differ distinctly from those found in other melanoma types. Most notably, the most common MAPK-activating mutations in melanoma, BRAF and NRAS, were completely absent in our cohort. Only a single tumor had an NRAS mutation, a rare type of RAS mutation in cutaneous melanoma. Three tumors had BRAF, BRAF, and BRAF substitutions – these mutations inactivate BRAF kinase activity but paradoxically activate MAPK signaling via activation of CRAF[35,36]. Intriguingly, the multifarious mutations in the RTK-Ras-MAPK/PI3K signaling cascade did not appear to be mutually exclusive, indicating that some of these alterations may cooperate in pathway activation (Fig 5A).

Figure 5

The mutational landscape of desmoplastic melanoma

A. Tiling plot of genetic alterations (rows) in each sample (columns). Mutations were considered homozygous if their MAFs were 1.5× the sample median. Numbers indicate the percent of samples harboring oncogenic alterations. Purple (*) tiles indicate unique hotspot mutations as discussed. GOF = gain-of-function. B. Pathways affected.

Other pathways with recurrent mutations (Fig 5) included the p53 and Rb pathways, mostly resulting from inactivation of CDKN2A, TP53, and RB1. Also, the SWI/SNF chromatin remodeling complex, a tumor suppressor in many malignancies[37], harbored inactivating mutations of ARID2 and ARID1A. TERT was activated in 90% of samples, mostly through promoter mutations[13,14] but also amplification.

The distinct landscape of genetic alterations in desmoplastic melanoma mirrors their unique clinical behavior. Some of the alterations found may have clinical implications. For example, there are small-molecule inhibitors directly targeting products of several oncogenes, such as MET, EGFR, ERBB2, IDH1, MAP2K1, PIK3CA, and CDK4. Furthermore, the exorbitantly high mutation burden found in desmoplastic melanomas makes them promising candidates for immune checkpoint blockade therapy[38]. Finally, the mutational patterns observed in desmoplastic melanoma indicate the ontogeny of these neoplasms. While desmoplastic melanomas typically present as deep-seated cutaneous neoplasms, the overwhelming UV-signature in most desmoplastic melanomas indicate that they arose from a cell in or near the epidermis that accumulated most of its mutation burden prior to dermal invasion.

Online Methods

Sample Collection and microdissection

The study protocol was approved by the UCSF's committee on Human Subjects, the Memorial Sloan Kettering institutional review board and human biospecimen utilization committee, and the Sydney Local Health District. Informed consent was obtained from all patients involved. Twenty fresh frozen desmoplastic melanomas and matching blood were acquired from Memorial Sloan Kettering (n=10) and The Melanoma Institute of Australia (n=10), comprising the discovery set. Forty-two archived FFPE desmoplastic melanomas and unrelated non-lesional “normals” were microdissected from the University of California San Francisco dermatopathology repository, comprising the validation cohort. There was no source of non-lesional tissue for 4/42 validation samples. All cases were initially diagnosed as desmoplastic melanoma and confirmed by an independent review. Tumors were assessed by H&E and adjacent unstained sections were manually microdissected by a pathologist to enrich for neoplastic cells. For DMs of the “mixed” subtype, best efforts were made to microdissect only the desmoplastic components. Clinical characteristics are summarized in table S1.

Sequencing, variant calling, and copy number analysis

DNA was extracted using Qiagen DNeasy kits and prepared for sequencing using the NuGen Ovation library preparation kit (p/n 0331-32) under manufacturer's protocols. Agilent SureSelect Exome V4 + UTR (p/n 5190-4638) or Nimblegen SeqCap EZ Libraries (p/n 06588786001) were used for capture of the exome or 293 genes (Supplementary Table 2) respectively. Multiplexed samples were sequenced on Illumina HiSeq 2500. Initial alignment was performed with Burrows-Wheeler Aligner (BWA), followed by INDEL realignment, deduplication, and recalibration by the Genome Analysis Tool-Kit (GATK). Point mutations were called using MuTect, and INDELs were called using Unified Genotyper and Somatic INDEL detector as compared to non-lesional “normals”. Somatic variants were annotated by Oncotator.

For the four samples in the validation cohort without non-lesional normal tissue (Supplementary Table 1), somatic mutations were inferred by searching for variants from the reference genome that have not been reported in 1000 Genomes or the ESP5400 NHLBI Exome Sequencing Project. To rule out private germline SNPs, mutations with mutant allele frequencies (MAFs) near 50% were removed. Because of stromal cell contamination, most somatic events clustered at MAFs less than 50%.

High resolution aCGH on Agilent 180K, 244K, or 1Mil feature arrays was performed on 44 samples in the cohort (Supplementary Table 1). Copy number was inferred from sequencing data in all 62 samples using the software package CNVkit (https://github.com/etal/cnvkit). Raw aCGH copy number data has been deposited in the Gene Expression Omnibus (GSE55150). Segmented copy number calls from aCGH and derived from CNVkit are included in table S5.

Melanoma subtype comparisons (related to Fig S1)

Genome-wide copy number profiles were compared between DMs and other melanoma subtypes[39] (Supplementary Figure 1). In order to do this, we needed to accurately call gains and losses in the presence of varying stromal cell contamination across samples. Towards this goal, we first used the scaled median absolute deviation (MAD) of the sequential order difference as implemented in the matrixStats package in R to estimate the sample-specific experimental variation. Next, we declared a segment to be gained or lost if the average log2 ratio was at least one sample MAD away from the median segmented values of the autosomal probes in that sample. Using this approach, the frequency of gains and losses are displayed for each melanoma subtype in the top panels of Fig S1.

We also explored whether other subtypes of melanoma had significantly different copy number alterations from desmoplastic melanoma. We tested association between copy number aberrations and melanoma subtype for each probe by performing 2-sided Fisher's exact test from trichotomous gain/loss/normal data. We next calculated q-values from these p-values with the Bioconductor q value package instituted in R. Significantly different regions of copy number alteration are displayed as indicated in the lower panels of Fig S1.

Calling TERT promoter mutations (related to Supplementary Figure 6)

TERT promoter mutations were discovered midway through the study. The region was tiled and sequenced for several of the latter samples in the study, but it was not directly sequenced for many of the initial samples. However, as the mutations reside near the first exon of TERT, there was low coverage over the TERT promoter even when it was not specifically tiled. For the discovery cohort, we could also supplement this coverage with whole genome sequencing reads. As there was a range of sequencing coverages and neoplastic cellularities, we calculated our power to call mutations over this locus for each tumor. Using a binomial test, we determined the likelihood that a true mutant sample with a given coverage and tumor cellularity could yield less than 2 mutant reads, thus yielding a false negative. We report the TERT promoter mutation status for all samples that could be called with 90% and 95% confidence (Fig S6).

Point Mutation and Copy Number Validation

The accuracy of mutation calling was accessed by resequencing a subset of mutations detected in the discovery phase to extremely high depth (>900× coverage). 336 mutations spanning 6 samples were selected for validation. In total, 329/336 (97.9%) of mutations were confirmed. These are noted in an extra column in Supplementary Table 3. In addition, 11 NFKBIE mutations were validated using Sanger sequencing as indicated in Fig S7. Finally, every sample exhibited specific mutational signatures that would not be detectable if the majority of mutations were artifacts.

Copy number alterations were initially inferred from array comparative genomic hybridization data. We additionally inferred CNAs from sequencing data using the CNVkit software suite[40]. The segmentation output from both approaches is included in table S5. Every single featured copy number alteration from Supplementary Figures 2-3 and Figure 3 was confirmed by both methods.

Nominating gene candidates for validation sequencing

Whole-genome or -exome sequencing was not feasible on the validation cohort of DMs due to the low tissue yields and fragmentation resulting from formalin fixation. The top candidate cancer genes (n=293 genes, listed in Supplementary Table 2) were nominated from the discovery set (n=20 samples) for targeted sequencing in a validation cohort (n = 42 samples). We included genes with >2 novel hotspot mutations in the discovery set. Any gene with a COSMIC hotspot mutation, even in a single sample, was also included. From the discovery set, we computed a loss of function burden for each gene and included the top 50 candidates for further validation sequencing. This calculation is described in the manuscript and in more detail below. Focally deleted and amplified genes were also included for validation sequencing. Finally, a curated list of known cancer genes was also included[41].

Calculating the loss of function burden

To identify genes with an elevated loss-of-function (LOF) mutational burden, we developed a computational method to compare the frequency and zygosity of loss-of-function mutations affecting each gene to that expected by chance. Our source code implementing this method is freely available at: https://github.com/etal/lofsigrank. A brief description of the computation follows.

Step 1: The overall non-synonymous:synonymous mutation ratio is calculated. In this study, that ratio was 2.13:1.

Step 2: The number of synonymous mutations striking each gene is multiplied by the overall non-synonymous:synonymous mutation ratio, as calculated in Step 1, and subtracted from the total number of mutations striking that gene. Genes for which this resulting number is positive are considered to be enriched for non-synonymous mutations and selected for further analysis.

Step 3: Within each gene, sample, and mutation type, non-synonymous mutation counts are multiplied by their normalized mutant allele frequencies (NMAF). The NMAF is calculated by dividing the mutant allele frequency of a mutation by the median mutant allele frequency of all mutations in a sample. Effectively, the NMAF approximates the clonality and zygosity of a given mutation. Thus, weighing mutations by their NMAF favored fully clonal mutations that had undergone loss of heterozygosity over subclonal or heterozygous mutations.

Step 4: For genes with multiple non-synonymous mutations in a single sample, the NMAFs of all mutations in that specific gene from that single sample are summed and capped at 2. The cap at 2 represents the maximum number of alleles that a single sample can lose. The rationale behind this approach is that multiple mutations striking the same gene in a single sample are likely to affect both alleles but can logically affect no more than both alleles. This approach avoids placing undue emphasis on large genes that accumulate high counts of incidental mutations. To illustrate the utility of this methodology, in this study, TTN, the largest gene in the genome, harbored over 100 mutations in a single sample – this correction credited that single sample with 2 hits to TTN (inactivating both alleles) rather 100+ hits.

Step 5: For each gene, sum the rescaled and capped mutation burdens across all samples, yielding a preliminary measure of the number of alleles impacted by mutations in a gene across all samples. This preliminary measure will be referred to as the number of hits.

Step 6: Determine the deleterious skew of non-synonymous mutations across all mutations in all samples. The frequency of truncating (nonsense, frame-shift, and splice-site) and missense (probably damaging, possibly damaging, or benign) mutations is calculated across all mutations in all samples. For the present study, the breakdown was 8.1% (truncating) and 29.6%/13.6%/16.7% for missense probably/possible/benign – the remainder were synonymous mutations.

Step 7: Determine the “truncating mutation” burden for each gene. The frequency of truncating mutations affecting a specific gene across all samples is calculated and divided by the genome-wide average. For example, if a gene has 5 truncating mutations out of 10 total mutations, then the truncating factor would be 0.5/0.081 = 6.17. This factor has a wide range of values.

Step 8: Determine the “missense-probably damaging” burden for each gene. The frequency of probably damaging mutations affecting a gene is calculated as follows: 1 + (missense-pd_gene) – (missense-pd_total). For example, if a gene has 5 missense probably damaging mutations out of 10 mutations total then this factor would be as follow: 1 + 0.5 – 0.296 = 1.204. This factor has a relatively small range of values.

Step 9: Calculate the loss of function burden as the product of the number of hits and mutation burdens calculated in steps 5 through 8: Multiply the normalized number of hits (step 5) by the truncating mutation burden (step 7) and the probably damaging mutation burden (step 8). As can be seen, the factors in steps 7 and 8 have the potential to magnify the value calculated in step 5 if the mutations affecting a gene were disproportionately damaging as compared to the genome average.

Step 10: Determine significance by simulation. Repeat steps 1-9 on sample/gene/mutation permuted data (recommended 1000 permutations). Calculate false discovery rates by comparing the observed loss of function burdens to permuted data.

NFKBIE PCR amplicon sequencing (related to Supplementary Figure 7)

Deep PCR amplicon sequencing was performed to validate a subset of mutations in NFKBIE (Supplementary Figure 7). 1ug of PCR amplicon was prepared for sequencing as described above and directly sequenced (without undergoing target capture). The resulting reads were aligned to the genome but were not deduplicated.

Discovery of NFKBIE mutant cases in non-desmoplastic melanoma

As part of ongoing efforts unrelated to this study, our group (BCB) is sequencing diverse pigmented lesions. Incidentally, we uncovered NFKBIE promoter mutations in 6 non-desmoplastic melanomas that were not explicitly part of this study. We've reported these mutations and the clinical features of those samples here (Fig 4, Supplementary Tables 1,6); however, the full body of work will be described in its entirety in future publications.

Inferring tumor cellularity

Tumor cellularity was inferred for each case and is listed in Supplementary Table 1. Tumor cellularity was estimated by doubling the median mutant allele frequency of somatic mutations from a given sample. This assumes that the median mutation allele frequency approximates a clonal, heterozygous mutation. Copy number alterations (CNAs) and subclonal mutations could skew this approximation. CNAs were infrequent in desmoplastic melanoma, so they would not be expected to impact our calculation. However, to more broadly confirm that CNAs and subclonal mutations had little effect on this calculation, we manually inspected the distribution of mutant allele frequencies in each sample. For each sample, we observed a predominant cluster of mutant allele frequencies centered around the median.

It is also possible to infer tumor cellularity from copy number alterations. We attempted to do this, but most samples were not adequate for analysis due to the scarcity of copy number alterations characteristic of desmoplastic melanoma. Notably, CNA derived estimates did agree well with our mutant allele frequency estimates for the subset of samples that were adequate for analysis.

RNA-Seq analysis of TCGA data (related to figure 4)

NFKBIE transcription was analyzed from RNA-Seq data of 60 melanoma TCGA (The Cancer Genome Atlas) samples. Only the short isoform of NFKBIE was expressed. The track in figure 4 was created by combining all 60 bam files into a single bam file and visualizing the cumulative read depth.

Immunohistochemistry and Western Blotting

Immunohistochemistry was performed using the following antibodies: p16 (Ventana, catalog number 9517), p53 (Dako, catalog number GA6166), Cyclin D1 (ThermoFisher, catalog number RM-914-S), CDK4 (Invitrogen, catalog number AHZ0202), EGFR (Ventana, catalog number 790-4347), MDM2 (Invitrogen, catalog number 182403), Rb (BD Biosciences, catalog number 554136), Met (clone c-28), YAP1 (Cell Signaling, catalog number 4912), phospho-ERK p44/42 (clone D13,14.4E).

Antibodies for Western Blots were as follows: NFKBIE (Santa Cruz, catalog number sc-7155); RELA/p65 (Cell Signaling, 8242); p50/p105 (Cell Signaling, p/n 3035); cREL (Cell Signaling, p/n 4727); PARP (Cell Signaling, 9542); HSP60 (Santa Cruz, sc-1722). NFKBIE siRNAs were from Dharmacon (ON-TARGETplus Human NFKBIE siRNA SMARTpool). Cell fractionation was performed using the Subcellular Protein FractionationKit for Cell Culture (Thermo Scientific 78840). The M257 and M375 cell lines were provided by Dr. Antoni Ribas.

Table S1. Sample Descriptions.

Table S2. Target Intervals. The 293 genes included in the validation cohort are listed in tab1. The second tab describes the actual genomic intervals tiled.

Table S3. Discovery Mutations.

Table S4. Validation Mutations.

Table S5. Copy number circular binary segmentation (CBS) calls. CBS calls[42] generated from the higher resolution primary copy number platform (Supplementary Table 1) are listed in the first tab. For many samples, copy number was also inferred using an orthogonal assay, and those copy number calls are listed in the second tab. Copy number data agreed very well between the two platforms.

Table S6. NFKBIE mutation details.