The molecular cartography of malignant and benign sebaceous tumours

Domain: biology
Published: 2025-12-19
Editorial outcome: major_revision

Abstract
Sebaceous tumours (STs) are rare skin appendage tumours and include benign sebaceous adenoma (SA) and sebaceoma (SM), malignant extra-ocular sebaceous carcinoma (SC-E) and peri-ocular sebaceous carcinoma (SC-O). Here, an extensive worldwide collection of 286 tumours is deeply characterised, revealing a propensity to develop in the context of a high tumour mutational burden (except in SC-O) which is most frequently associated with mismatch repair deficiency (dMMR), followed by UV-induced damage, POLE/POLD1 mutations, and AID/APOBEC activation signatures. Biallelic TP53 inactivation with concomitant ZNF750 and/or RB1 mutation is seen in SC-E/SC-O. Amplification of 8q (including MYC) is related to SC-O, while amplification of 1q21.3 (including HRNR) and chromosome 20 are shared by SC-O and SC-E, as is deletion of 13q14.3 (where RB1 resides). The most frequently mutated gene is NOTCH1. Extensive fusion gene, expression and molecular cluster analyses provide a molecular portrait of this rare and enigmatic tumour type.

Introduction
Sebaceous tumours (STs) are rare cutaneous and peri-ocular neoplasms of sebaceous gland origin, including benign sebaceous adenoma (SA) and sebaceoma (SM), and malignant sebaceous carcinoma (SC). SC is anatomically divided into peri-ocular (SC-O) for tumours located on the eyelid (conjunctiva and/or skin aspect), and extra-ocular (SC-E) for tumours elsewhere on the skin. Clinically, SA and SM typically present as solitary and asymptomatic slow-growing yellow nodules (occasionally ulcerated) on adult seborrhoeic skin areas and are treated by complete excision1,2. SC commonly appears as a solitary nodule in adults of 60–79 years3,4, with SC-O more common in people of Asian descent (40–60% of eyelid tumours). SC-O and SC-E are associated with significant morbidity and mortality, showing recurrence in 11.8% and 14% of patients, metastasis in 12.5% and 1.8% of cases, and a disease-specific mortality of ∼25.6% and ∼18%, respectively3. The recommended treatment is complete excision (including orbital exenteration), with radiotherapy used for inoperable cases or as an adjunct to surgery3.

STs can develop sporadically, as mismatch repair (MMR) proficient (pMMR) or deficient (dMMR) tumours, and may occur in association with Lynch Syndrome (LS), which results from constitutional (germline) pathogenic/disruptive variants in one of four MMR genes; MSH2, MSH6, MLH1, or PMS2, with LS-STs deficient for MMR (dMMR) also exhibiting microsatellite instability (MSI)5. LS patients have an increased lifetime risk of developing colorectal (CRC) and endometrial cancer, as well as a wider range of other primary cancers, including STs and keratoacanthoma. Muir-Torre Syndrome (MTS) is a subtype of LS that is characterised by at least one visceral LS-associated cancer and one cutaneous sebaceous tumour (except SC-O)6.

We perform a large-scale analysis of the sebaceous spectrum integrating whole-exome and transcriptome sequencing data7–13 of an extensive globally ascertained cohort including SAs, SMs and SCs totalling 286 cases to provide a comprehensive genomic portrait of these tumours and the molecular events that drive their development.

Results
Sample ascertainment and clinico-histopathological features
Samples from 222 patients (286 tumours, of which 197 primary ST had matched/paired normal DNA) were ascertained from 11 institutions across six countries (Belgium, Canada, Germany, Greece, Japan, United Kingdom). To explore the genetic origin of the patients, their genotypes were projected against the 1000 genomes dataset (Supplementary Fig. 1 in Supplementary Information), revealing largely European and Asian descent as expected. Formalin-fixed paraffin-embedded (FFPE) samples were comprised of four ST subtypes, specifically SA (n=102), SM (n=92), SC-E (n=49 primaries and n=1 metastasis), and SC-O (n=26 primaries, n=2 recurrences and n=7 metastases); the latter being underrepresented due to the rarity of this ST subtype, which might limit the power of certain analyses. Additional related lesions were included in the study for comparative purposes, specifically sebaceous hyperplasia (SH, n=6) and a keratoacanthoma with sebaceous differentiation (KA, n=1). For the four ST subtypes (SA, SM, SC-E and SC-O), the median age was 73, 69, 76 and 75 years at diagnosis, with a sex ratio M/F of 2.5/1, 1.04/1, 1.9/1, and 1.2/1, respectively. Most tumours were located on the head and neck (65.8%), followed by the trunk (16.4%), the peri-ocular area (12.6%, of which 9.7% were SC-O), the limbs (3%) and the genital area (1.9%), while one case was from an unspecified site (0.3%). A summary of the clinico-histopathological features of the four ST subtypes, as well as an overview of the project is shown in Fig. 1A–E (Supplementary Data 1 and Supplementary Discussion in Supplementary Information).

Immunohistochemistry (IHC) for MMR proteins (MSH2, MSH6, MLH1 and PMS2) was performed on 235 tumours (on at least one tumour per patient). The results of the MMR IHC were used together with indel (insertion/deletion) rate, dMMR mutational signatures, and somatic and germline variant calls to classify the MMR status of tumours (pMMR and dMMR; see “Methods”, Supplementary Data 1 and Supplementary Discussion in Supplementary Information).

Immunosuppression is associated with an increased incidence of STs, and acquired dMMR can be recognised in some of these tumours14. Three out of the six tumours from immunocompromised patients were dMMR in our cohort. Other neoplasms occurred in some patients and were found in individuals who had presented with any of the four ST subtypes (n=103; Fig. 1B, Supplementary Data 1). ST development can be the first clinical sign of LS seen in from 9% and possibly up to 52% of the 23 LS patients in our cohort. Critically, this highlights that the diagnosis of an ST could be used to initiate appropriate cancer screening or potentially prophylactic treatment strategies15. To evaluate this possibility further, the clinical and histopathological attributes associated with LS were examined in our cohort using univariable and multivariable logistic regression. Of 6 clinical attributes, developing multiple STs, having a tumour outside the head and neck, and having another non-ST LS-related neoplasm were all significantly associated with LS after multivariable testing and following false-discovery correction (multivariable odds ratios of 22, 14, 8.8, respectively; Fig. 1E and Supplementary Discussion in Supplementary Information). Ten histopathological features were also tested. No epidermal connection and delimitation of the ST by a collarette appeared to be statistically significantly related to LS on multivariable testing and after false-discovery correction (multivariable odds ratio 5.6 and 4.8, respectively; Fig. 1E and Supplementary Discussion in Supplementary Information).

Follow-up was obtained for 179 patients (mean: 39 months, median: 26 months; range: 1–144 months), and from the clinical data it was observed that SC-O patients had a worse prognosis than SC-E patients; recurrence and metastasis occurred in 7.7% and 19.2% of the SC-O patients after an average time of 8.5 and 48 months compared with 2.2% and 8.7% of the SC-E patients after an average time of 21 and 16.2 months, respectively (Supplementary Data 1). Of interest, metastasis occurred to lymph nodes for SC-E, while lymph node, parotid gland, lung and brain metastasis was observed for SC-O. Tumour-related death occurred in 7.7% and 4.3% of the SC-O and SC-E patients, respectively. A summary of the individual sample details is provided in Supplementary Data 1.

Somatic landscape of sebaceous tumours
Whole-exome sequencing (WES) was performed to profile somatic mutations, somatic copy number alterations (SCNAs) and mutational signatures. The tumour mutational burden (TMB; including single nucleotide variants (SNVs), multi-nucleotide variants (MNVs) and insertions/deletions (indels)) varied markedly across the cohort. Using mutation calls from the 197 primary STs where a germline matched normal sample had also been sequenced, we found a statistically significantly higher TMB for dMMR STs compared to all other tumours in The Cancer Genome Atlas (TCGA) dataset, and dMMR STs also had a much higher TMB compared to pMMR STs (linear model with pairwise comparisons using Tukey’s method; Fig. 2A–C and Supplementary Data 2). dMMR STs also showed a higher TMB than a range of skin tumours sequenced by MSK-impact (linear model with pairwise comparisons using Tukey’s method; Supplementary Data 2). Interestingly, indel rate (Area Under the Curve; AUC 0.999) was a better predictor of dMMR status compared with TMB (AUC 0.887), as there were pMMR tumours with a low indel rate and high TMB (Supplementary Fig. 2A in Supplementary Information). Importantly, for this specific analysis, we defined dMMR cases as those having lost one of the MMR proteins by immunohistochemistry and also carrying biallelic inactivating events in an MMR gene. We did not use the indel or SNV rate (see “Methods”). Using a threshold of 0.1 on the predicted probability of dMMR computed with our model, a rate of >0.345 indels/Mb will give 97% sensitivity and 99% specificity to predict dMMR. However, this threshold is specific to our pipelines, as the sensitivity and accuracy of mutation calling is an inherent property of each variant calling pipeline. Amongst the different ST subtypes, SC-Es had the highest TMB regardless of MMR status, while SC-Os (which were all pMMR) had a TMB comparable to that of benign pMMR tumours (Supplementary Data 3). Interestingly, pauci-mutated samples (with <5 mutations/Mb) were found across all four ST subtypes, usually related to pMMR status, however, in SC-O, pauci-mutation was the predominant phenotype (Supplementary Fig. 2B in Supplementary Information and Supplementary Data 3).

The complete list of somatic variants is provided in Supplementary Data 4 (with cross-validation using RNAseq data provided in Supplementary Fig. 4A in Supplementary Information and see “Methods”). Considering only the paired primary tumours (n=70 SA, n=70 SM, n=38 SC-E and n=19 SC-O; n=197 in total), the top 20 recurrently mutated genes mainly consisted of known cancer-associated genes (Supplementary Fig. 3A, B in Supplementary Information). Driver gene analysis of the somatic variants across all four ST subtypes (using dNdScv; see “Methods”) identified 25 statistically significantly mutated genes (Fig. 2A, Supplementary Data 5 and Supplementary Fig. 3C in Supplementary Information includes the unpaired tumours). Additional driver genes were identified when considering only specific ST subtypes and/or MMR status, such as PRRT2 in dMMR tumours, DRD5 in dMMR SC-Es and SPANXB1 in dMMR-MSH2/MSH6-negative SC-Es (Fig. 2A, Supplementary Fig. 3C in Supplementary Information and Supplementary Data 5A). To better understand the role of these driver genes in STs, the variant types they carry (Supplementary Fig. 4B in Supplementary Information), their presence in SCNAs (Supplementary Fig. 4C in Supplementary Information) and their predicted mechanism of contribution to tumorigenesis (oncogene versus tumour suppressor gene (TSG); Supplementary Fig. 4D in Supplementary Information), and their variant allele frequency (VAF; Supplementary Fig. 4E in Supplementary Information) were considered. HRAS, TMEM86A and SGCB were characterised by missense mutation only, suggesting they are acting as oncogenes; while most driver genes show a range of variant types (including indels and truncating mutations), more suggestive of a TSG profile16. Similarly, several established oncogenes were amplified (such as KRAS), while loss and copy neutral loss of heterozygosity (cnLOH) were identified at tumour suppressor gene loci, such as MSH2 and RB1. In addition, these 28 driver genes were cross-referenced with pan-cancer analyses17, to better understand their impact and activities in other cancers (Supplementary Fig. 4F in Supplementary Information); 15 out of the 28 driver genes (53.6%) identified by dNdScv were shared with genes in this pan-cancer driver gene resource. Additionally, reduced hypothesis testing was also performed with 312 previously defined cancer driver genes from Bailey et al.17, allowing us to identify further potential driver genes (Supplementary Fig. 5 in Supplementary Information and Supplementary Data 5B). Lastly, using another tool (OncoDriveFML), further candidate driver genes were called (see GitHub repository), while 22 out of the 25 drivers genes (88%) found by dNdScv were also found with this method.

The transmembrane receptor NOTCH1 was the most recurrently mutated driver gene (n=112/197), with many tumours having multiple hits (n=67/112). The mutations predominantly occurred in the extracellular epidermal growth factor (EGF) domain (n=87/112; 70/87 of these tumours were dMMR), with a hotspot at p.A465T identified (n=10; COSM1217602; Supplementary Fig. 6A in Supplementary Information). NOTCH1 can function as both a tumour suppressor gene and oncogene, depending on the cellular context18, however, the mutations observed in STs were primarily inactivating and this suggests that NOTCH1 predominantly functions as a TSG. The next most recurrently mutated driver gene was the transcription factor RAS responsive element binding protein-1 (RREB1; n=86/197; Supplementary Fig. 6B in Supplementary Information), with many tumours having multiple mutations (n=47/86). Tumours with RREB1 mutations predominantly (but not exclusively) also had NOTCH1 mutations (n = 71/86). The HRAS oncogene was the third most recurrently mutated driver gene (n=81/197), with predominantly missense mutations in the well-established hotspots of exon 2, p.G12C/D/R/S (n = 19/81) and p.G13D/R/S/V (n=17/81), and exon 3, p.Q61H/K/R (n=10/81), as well as other less well known variants that are reported in COSMIC, such as p.K117N (n=12/81), p.A18V (n=11/81) and p.A146T (n=5/81; Supplementary Fig. 6C in Supplementary Information and Supplementary Data 4). Interestingly, when considering the ST subtypes individually, HRAS was only identified as a statistically significant driver gene in the benign tumours. Notably, hotspot HRAS mutations were observed in malignant subtypes, but these mutations occurred at a much lower mutation frequency in these smaller cohorts (Supplementary Data 5A). In addition, specific variants of HRAS were found to associate with MMR status, as pMMR tumours had p.G13R/V or p.Q61K/R variants, whilst dMMR tumours specifically harboured p.G12S, p.A18V, p.R68Q/W and p.A146P/T/V variants (with the p.K117N/R variants being associated with either MMR status). Another gene found to associate with MMR status was the BAP1 complex regulator methyl-CpG-binding domain protein 6 (MBD6), in which mutations (predominantly frameshift insertions or deletions, with hotspot mutations at p.A783Sfs*10 (n=24), p.G782Efs*13 (n=17), p.Q694Sfs*5 (n=10) and p.Q734Nfs*18 (n=6)) were only found in dMMR tumours (n=62/197, except for one pMMR SC-E). The fourth most recurrently mutated driver gene was TP53 (n=79/197), with missense mutations predominating (n=36/79). The presence of mutations in the top four most mutated driver genes allowed the identification of five molecular subtypes, specifically those with NOTCH1 mutations in the presence of RREB1 mutations (molecular subtype A), NOTCH1 mutations in the absence of RREB1 mutations (molecular subtype B), RREB1 mutations in the absence of NOTCH1 mutations (molecular subtype C), HRAS mutations in the absence of NOTCH1 or RREB1 mutations (molecular subtype D), and TP53 mutations in the absence of NOTCH1 and RREB1 mutations (molecular subtype E); highlighting a shared mutational profile amongst the four ST subtypes despite their histological differences (Fig. 2A). Whilst molecular subtypes A, D and E could be identified in SC-O, they also carried mutation of TP53 (n=18/19) and the transcription factor ZNF750 (n=13/19; TP53 and ZNF750, n=12/19; Fig. 2A), which was also seen in the metastatic SC-Os (Supplementary Fig. 3B in Supplementary Information). This is consistent with a recent exome study of five eyelid SCs19. A proportion of SC-Es also had mutations in both TP53 and ZNF750 (n=10/38). Conversely, the dopamine receptor 5 (DRD5) gene was mainly identified as a driver gene in dMMR SC-Es associated with both NOTCH1-RREB1 mutation (mutated in 9/38; Supplementary Data 5A). Finally, both SC-E and SC-O showed recurrent mutations in the TSG, RB1 (n=8/38 and n = 6/19, respectively), in contrast to a previous report in which RB1 mutations were only found in SC-O (not SC-E)10. Of note, both ZNF750 and RB1 mutations were largely associated with TP53 mutation in SC (in 85% and 86% of cases, respectively, of SC-E and SC-O cases combined (Fig. 2A, B). In addition, concomitant TP53/RB1 mutations were identified in metastatic SCs. When stratifying by ST subtypes, and using DISCOVER20, three mutually exclusive interactions were identified (ZNF750 remained mutually exclusive from RREB1, MIDEAS and BICRA (q-value=0.038 for each)), which were also seen amongst the 17 mutually exclusive interactions observed in the unstratified analysis. Unstratified DISCOVER analysis (see “Methods”) also found that NOTCH1 (q-value=0.015), RREB1 (q-value=0.008) and HRAS (q-value=0.001) mutations were mutually exclusively from TP53 mutations, and NOTCH1 (q-value=0.010) mutations were mutually exclusively from ZNF750 mutations, confirming the separation of the STs into five molecular subtypes, especially as molecular subtype E has additional driver gene mutual exclusivity (Fig. 2B and Supplementary Data 6).

Integrative cluster analysis combining the somatic mutation and transcriptional profiles (see GitHub repository) of the 170 STs for which both exome/DNA and transcriptome/RNA analyses had been performed (see “Methods”), identified 9 integrative clusters (IC; Fig. 3, Supplementary Fig. 7 in Supplementary Information and Supplementary Data 7), each characterised by distinct clinico-pathological features, indel/SNV rates, somatically mutated genes, gene expression/pathway enrichment, and tumour microenvironment (TME) features. Interestingly, these individual ICs shared many features with the individual molecular subtypes. For example, IC1–IC3 and IC6 overlapped with molecular subtype A (standardised residuals 2:4 each), which is characterised by a high proportion of mutations in driver genes such as NOTCH1 and RREB1 (with or without BICRA and MIDEAS mutations) and was found mostly in benign dMMR tumours. They also contained a high proportion of ASXL1 mutations in addition to NOTCH2 in IC3 and IC6, and MBD6 in IC6. Interestingly, these three ICs, driven either by mutational burden (IC1 and IC6) or indel burden (IC3), showed different pathway enrichments and a different tumour microenvironment, with IC1 having up-regulation of G2M checkpoint and E2F targets, and containing a high proportion of CD8+ T cells and fibroblasts, IC3 being down-regulated of most hallmark pathways and showing a high proportion of CD4+ T cells, and IC6 being up-regulated in the MYC pathway; suggesting possible differences in the tumour behaviour or cellular makeup in each molecular subtype. IC7 showed similarities with NOTCH1 mutated molecular subtype B (standardised residuals 2:4). Displaying a lower mutational burden profile, this IC (IC7), mostly encountered in pMMR and dMMR SAs, is mutated for NOTCH1 and HRAS in 50% of cases, showed a high proportion of endothelial cells and fibroblasts, and is down-regulated for interferon pathways. Interestingly, IC8 differed from molecular subtype A (standardised residuals -4: -2) and overlapped with molecular subtype D (standardised residuals >4). Mostly composed of pMMR benign STs, IC8 is characterised by a low mutational burden with mutations in HRAS, ZNF750 and TP53 occurring in 40 to 50% of cases, the contribution of UV related mutational signatures, and down-regulation of many pathways. Finally, molecular subtype E was linked to IC4 and IC9 (standardised residuals 2:4 and >4, respectively), both showing a high proportion of mutation of TP53 and larger copy number alterations (at a chromosomal arm level) compared to the other ICs. IC4 exhibited additional ZNF750 and/or NOTCH1 mutations in 40% of cases and was mostly composed of pMMR SCs (SC-Es predominantly), included a high population of B and T cells, showed UV and APOBEC related mutational signatures, and up-regulation of the interferon response. In contrast, IC9 cases had an even higher proportion of mutations in ZNF750 with some cases also mutated for RB1. These cases were exclusively pauci-mutated (pMMR) SC-Os and showed up-regulated of E2F and G2M checkpoint pathways, and a lower immune component.

Mutational signatures of sebaceous tumours
Mutational signatures are characteristic patterns of somatic single-base substitutions (SBS), doublet-base substitutions (DBS) and small indels (ID) that reflect underlying mutational processes21. As expected, in the dMMR cases, the high TMB was attributed to dMMR-associated mutational signatures, specifically SBS15, SBS20, SBS21, SBS26 and ID7 (Fig. 4A, B, Supplementary Fig. 8 in Supplementary Information and Supplementary Data 8), consistent with previous reports10,22. Interestingly, most tumours with ID7 signatures also had ID2 signatures; ID2 is associated with slippage during DNA replication of the template DNA strand and tends to be highly elevated in dMMR cancer samples21. Signatures associated with ultraviolet (UV) light exposure, specifically SBS7a, SBS7b, SBS7d and DBS123,24, were also associated with a high TMB and were identified in a proportion of samples from all four ST subtypes (DBS1 being identified in SC-E), validating previous observations10,22. However, SBS7 was most notably found in the SC-Es (n=16/38, of which 6 also showed a DBS1 mutational signature; all were from sun-exposed sites, and n=14/16 were pMMR). Importantly, one case (PD42510a) of the 19 SC-Os harboured UV signatures, and while all other SC-Os were from the eyelid or eye, this sample was from the canthus (Supplementary Fig. 9 in Supplementary Information). Of note, signatures SBS2 and SBS13, which are attributed to the activity of the AID/APOBEC family of cytidine deaminases, were found in some pMMR SC-Es (n=5/38; all of which were from sun-exposed sites and pMMR; Fig. 4A, B).

Interestingly, one sample with dMMR-associated signatures and a high SNV/indel rate (PD45524a; Fig. 4A) was also found to have signatures SBS10a and 10b, which are associated with POLE exonuclease proofreading defects. In CRC, POLE mutation results in a large number of somatic mutations (>100 mutations/Mb) and are described as ultramutated25 as compared to the dMMR-associated hypermutated samples21. Importantly, this sample had somatic mutations in the exonuclease domain (ED) of POLE (n=2) and POLD1 (n=1), as well as in MSH2 (n=1) and MSH6 (n=3) and 273 mutations/Mb. Four other samples had one somatic mutation each in the ED of POLE (PD43395a, PD43745a, PD43910d, PD47239a), however, SBS10a/b mutational signatures were only identified in one of them (PD43745a). Similarly, signature SBS20, which is associated with dMMR in the context of concurrent POLD1 mutation26, was identified in seven tumours, however, in four tumours a POLD1 mutation could not be found (PD43911a, PD47202d, PD42500a, PD43869a). STs associated with somatic POLD1 and POLE mutations and/or their related mutational signatures have never been reported before. One sample was found to have signature SBS36, which is correlated with defective DNA base excision repair due to disruption of MUTYH (germline or somatic)27, and a homozygous MUTYH pathogenic allele was identified in this case (germline; PD42485c). Interestingly, signature SBS32, which is linked to prior treatment with the immunosuppressive drug azathioprine28, was attributed to a SC-E from a renal transplant patient (PD42480a).

Copy number landscape of sebaceous tumours
Analysis of somatic copy number alterations (SCNAs) revealed distinct profiles between the four ST subtypes (Fig. 5). The benign tumours had relatively few chromosomal aberrations (Fig. 5A, B), whereas the malignant tumours were replete with large chromosomal gains and losses (Fig. 5C, D). There was overlap in the SCNA profiles of the malignant tumour subtypes, with both showing recurrent copy number (CN) gains of chromosomes 1, 8q and 20 and CN loss of chromosomes 4, 8p and 13q, however, the frequency of these CN aberrations was higher in SC-O than SC-E (Fig. 5C, D).

When taking all four ST subtypes together and separating them based on MMR status, SCNAs predominantly occurred in the pMMR tumours (Fig. 5E), as the dMMR tumours with loss of MSH2/MSH6 or MLH1/PMS2 expression by immunohistochemistry showed relatively few SCNAs (Fig. 5F, G). This suggests two different types of genomic instability occur with little overlap, specifically chromosomal instability in pMMR tumours and microsatellite instability (MSI) in dMMR tumours. Of note, regions of recurrent cnLOH differed based on the tumours MMR status. For example, recurrent chromosome 17 cnLOH was observed in the pMMR cases (Fig. 5H), with TP53 cnLOH seen in SC-Os (n=5/15), SC-Es (n=6/36) and SAs (n=3/68), and ZNF750 cnLOH observed in SC-Os (n=8/15), SC-Es (n=6/36), SMs (n=11/68) and SAs (n=2/68) (Supplementary Data 9). TP53 and ZNF750 cnLOH resulted in biallelic inactivation in 100% and 81.5% of cases, respectively, as these tumours also had a somatic mutation in these genes. Importantly, TP53 inactivation by cnLOH has been reported in several cancer types29–34; it has never been reported in STs. Recurrent MSH2 and MSH6 cnLOH was identified in 38/187 tumours (Fig. 5I), of which 89.5% were dMMR with loss of expression of MSH2/MSH6 by immunohistochemistry. Of these tumours, 66% were from LS patients with a constitutional (germline) disruptive/pathogenic MSH2 variant (Supplementary Data 9). Similarly, MLH1 cnLOH was observed in 11 tumours (Fig. 5J), out of the 21 dMMR tumours with a loss of MLH1 expression, with one tumour coming from an LS patient. Of note, two pMMR SCs (one SC-E and one SC-O) with no loss of MMR expression by immunohistochemistry exhibited MLH1 cnLOH. PMS2 cnLOH was identified in one pMMR SC-O. MMR inactivation through MLH1 and MSH2 cnLOH have been reported in MSI CRC cell lines and/or tumours35–37. Lastly, chromosome 9q cnLOH, including CDKN2A (n=8/187) and/or NOTCH1 (n=20/187; mostly mutated for NOTCH1), was identified predominantly in pMMR tumours, and to a lesser extent in dMMR tumours (Fig. 5H and Supplementary Data 9). Interestingly, 9q cnLOH has been previously reported in normal skin38.

To identify candidate driver genes in regions of CN gain or loss, both focal regions (≤ half a chromosome arm; Supplementary Data 10) and broad regions (> half a chromosome arm; Supplementary Data 11) were considered. Focal peaks were defined as significant when reaching the 4 criteria: (i) GISTIC2 residual q-value<0.1, (ii) peak size ≥100,000 bp, (iii) concordance with ASCAT CN calls ≥0.75, (iv) portion of peak not overlapping with Genome In A Bottle difficult regions by >40% (Supplementary Data 10); while broad peaks were determined as significant when reaching the 3 criteria: (i) GISTIC2 residual q-value<0.1, (ii) concordance with ASCAT CN calls ≥0.75, (iii) gain with log ratio >0.25 or loss with log ratio <−0.25 (Supplementary Data 11 and see “Methods”). There were no significant focal or broad SCNAs shared between the benign and malignant tumours. SC-Os showed the greatest number of focal SCNAs (Fig. 6). Oncogenes in these regions included GNAS, NFATC2, PTK6 and SALL4 (n=14/15 SC-Os), and LYN and PLAG1 (n=14/15 SC-Os). Deletions were also observed in the TSGs EZH2, TRIM24, CNTNAP2 and KMT2C (n=5/15 SC-Os). Of note, the focal SCNAs in the SC-Es only included known TSGs, specifically FOXO1 and RB1 (n=8/36), and ARHGEF10 (n=9/36). Interestingly, a significant focal amplification was seen at 1q21.3 in SCs (n=17/36 SC-Es and n=10/15 SC-Os), which contains HRNR (Hornerin, a S100 family protein; Fig. 6 and Supplementary Data 10). Similarly, a significant focal deletion at 13q14.13 was observed in 11/15 SC-Os, including no known TSGs (Fig. 6 and Supplementary Data 10).

At a chromosome arm level, SC-Os showed the greatest number of significant recurrent SCNAs, including chromosomes 1q, 8q, 19p and 20 amplifications, as well as chromosome 4 and 13q deletions (Supplementary Fig. 10A–D in Supplementary Information), with the 13q deletions encompassing RB1 and BRCA2 (n=8/15; one of them being somatically mutated for RB1; Supplementary Fig. 10E in Supplementary Information). In contrast, SC-Es showed significant chromosomes 1 and 20 amplification, and significant 17p deletions, which encompassed TP53 (n=7/36 tumours, of which 6 also had a somatic TP53 mutations; Supplementary Data 11).

Germline alleles in sebaceous tumour patients
Putative pathogenic constitutional (germline) variants (SNVs and indels) were inspected in each of the four ST subtypes (n=154 patients/197 samples; Supplementary Data 12). To focus on clinically relevant genes, only genes known to be associated with cancer-predisposition were examined (see “Methods”) and the variants were filtered to consider only disruptive variants, specifically, frameshift, nonsense and splice site variants (as annotated by Variant Effect Predictor (VEP)39) and/or alleles reported as pathogenic/likely pathogenic (as defined by ClinVar40). Using these criteria, the most recurrently mutated gene was MSH2 (n=12/154; 8% patients) with ClinVar clinical significance defined pathogenic/likely-pathogenic variants in 11 patients (Fig. 7). The remaining patient carrier of an MSH2 germline variant, which has been reported as non-pathogenic/likely pathogenic in ClinVar, developed multiple STs as well as Lynch-associated internal cancers (PD43912), suggesting the variant (p.Ser153ProfsTer21) in this individual may be disruptive and thus disease causing. Interestingly, two patients carried ≥2 variants in MMR genes; patients with variants in more than one MMR gene are termed digenic/double heterozygous LS, and while it is not yet clear if this results in a more severe form of LS, there are clear implications for genetic counselling41,42. Lynch-like syndrome (LLS) patients develop LS-associated neoplasms (in particular CRC) but do not carry constitutional variants in MMR genes, or somatic BRAF mutations or MLH1 hypermethylation in CRCs43 they develop. A proportion of these cases are caused by biallelic MUTYH inactivation44,45. In our ST cohort, four patients carried a disruptive/pathogenic constitutional variant in MUTYH (Fig. 7). Three of these tumours were pMMR; two of the patients developed colonic neoplasia (adenoma or CRC), consistent with the LLS phenotype, one of them (PD42485b) having a homozygous genotype.

Lastly, other ClinVar disruptive/(likely) pathogenic variants were identified in several cancer predisposition genes, including BRCA2 (n=6 patients), MSH6 (n=5 patients), MLH1 (n=4 patients) and TP53 (n=2 patient). Of note, no PMS2 constitutional variants were found.

Identification of fusion genes in sebaceous tumours
In addition to WES, exome capture RNA-seq was performed on tumours for the purpose of the IC analysis (see above) and for identifying fusion genes in the transcriptome and viruses and other pathogens/organisms from unmapped reads (Supplementary Discussion in Supplementary Information, Supplementary Data 13 and see “Methods”). 105 fusions were identified across all samples, including 9 that were recurrent between tumours from different patients or were found in distinct (unrelated) tumours from the same patient. Of these, 4 were inframe fusion events (Supplementary Data 14). Notably, PAK2 fusions (PAK2::DLG1, PAK2::LRIG1 and DLG1::PAK2) were identified (Supplementary Fig. 11A in Supplementary Information and Supplementary Data 14) with the reciprocal event DLG1::PAK2 observed in the case with a PAK2::DLG1 fusion. We did not identify a reciprocal event in the tumour with the PAK2::LRIG1 fusion, which might suggest it is expressed below the level of detection or it is not generated. Importantly, our analysis suggests that PAK2 fusions are associated with benign tumours and with sebaceous differentiation, as was recently reported in apocrine poroma46, with fusions found in two SMs and one SA, without specific morphological features, in our cohort (Supplementary Data 14). Previously unreported fusions (Supplementary Discussion in Supplementary Information) included RCOR1::GRHL2, an interchromosomal fusion identified in two sebaceomas from the same patient, both showing the same breakpoint (Supplementary Fig. 11B in Supplementary Information), suggesting a shared, possibly post-zygotic, origin with this specific patient, noted as having multiple recurrent sebaceous cysts on his back.

Predicting the immune landscape of sebaceous tumours
Factors that restrain ST tumour growth and malignant progression likely include contributions from the immune system, as seen in melanoma and other cancers. dMMR cancers carry numerous tumour neoantigens, which are critical targets of the host anti-tumour response and correlate with the efficacy of immunotherapy. Thus, tumour neoantigens were predicting across the cohort (see “Methods”). Hypothetical predicted neoantigens were analysed at the gene and variant level. Across the cohort, the top three genes harbouring mutations predicted to produce neoantigens (combined class-I and class-II) were the driver genes HRAS, RREB1 and NOTCH1, with neoantigens predicted in 28, 20 and 19 samples, respectively (Supplementary Fig. 12A in Supplementary Information and Supplementary Data 15). At the variant level (Supplementary Fig. 12B in Supplementary Information and Supplementary Data 15), HRAS showed recurrent predicted neoantigens (combined class-I and class-II) at p.G12S (n=7), p.K117N (n=4), p.A18V (n=3), p.Q61K (n=3), p.Q61R (n=2) and p.G13R (n=2; Supplementary Fig. 12B in Supplementary Information and Supplementary Data 15). ACVR2A p.K437Rfs*5 was the most recurrent variant (n=11), while NOTCH1 showed 19 distinct predicted neoantigen variants. The predicted neoantigen RREB1 p.P1064Rfs*19 was found in three samples. Other driver genes with putative neoantigens included ZNF750, MBD6, NOTCH2, PTPN14 and TP53 (Supplementary Data 15).

Taking all four ST subtypes together or separately, predicted class-I and class-II neoantigen load, as well as hypothetical predicted immunogenicity of the neoantigens, were significantly higher in dMMR tumours compared with pMMR tumours (Supplementary Fig. 13A–H in Supplementary Information). However, no differences between dMMR and pMMR tumours was observed for the SC-E cohort, apart from a minor increase in class I neoantigenicity. Interestingly, the SC-E cohort (regardless of MMR status) showed similar neoantigen counts and neoantigenicity as dMMR SA and SM. Looking at the subgroups identified by the Integrative Clustering analysis, IC3 and IC5 had the highest neoantigen tumour immunogenicity class-I, while IC4, IC8 and IC9 were among the lowest (Fig. 3D). In comparison, IC1 to IC7 were in the highest groups of neoantigen tumour immunogenicity class-II and combined, while IC9 was among the lowest part.

Methods
Cases election and material extraction
Ethical approval for the use of these samples and associated data was obtained by a local committee at the institution of origin and via Research Governance at the Wellcome Sanger Institute. This study is part of a larger study that has been approved by the National Health Service Health Research Authority (Research Ethics Committee reference 21/PR/1024, IRAS project ID 304621).

The samples (SA (n=102 primary tumours), SM (n=92 primary tumours), SC-E (n=49 primary tumours and n=1 metastasis), and SC-O (n=26 primary tumours, n=2 recurrences and n=7 metastases), SH (n=6) and KA with sebaceous differentiation (n=1)) were identified from eleven clinical centres across six countries. Representative haematoxylin & eosin (HE)-stained sections were independently reviewed by two specialist dermatopathologists to confirm diagnoses and identify areas for sampling (tumour+/−normal tissue). Immunohistochemistry stains were also reviewed when possible or were performed to aid diagnosis. Clinical data and follow-up were obtained from patient records in keeping with ethical approvals.

All tumour and normal tissue samples were obtained from formalin-fixed paraffin-embedded (FFPE) tissue blocks by either using a tissue microarray (TMA) coring needle or macrodissecting unstained 10-micron thick tissue sections attached to glass slides to collect specifically identified tumour and normal material following a detailed assessment of the H&E slide of each case. Genomic DNA and RNA was extracted from the tumour samples (with genomic DNA only extracted from normal samples) using the AllPrep DNA/RNA FFPE Kit (Qiagen), according to the manufacturer’s instructions.

Sequencing and data processing
Whole exome. Sequencing libraries were prepared from the FFPE-extracted DNA using a NEBNext Ultra II DNA Library Prep Kit (New England Biolabs), according to the manufacturer’s instructions. Unique dual index tags were applied, and the samples were amplified by PCR using the KAPA HiFi Kit (KAPA Biosystems) for a minimum of eight cycles. The libraries were quantified using Accuclear dsDNA Quantitation Kit (Biotium), pooled (8-plex) in an equimolar fashion and hybridised overnight with the SureSelect Human All Exon V5 baits (Agilent). The multiplexed samples were paired-end sequenced using the HiSeq and NovaSeq platforms (Illumina) to generate 101bp reads. Sequencing reads were aligned to the GRCh38 reference genome, using BWA-MEM50, and PCR duplicates were marked using Biobambam2 bammarkduplicates2 (v2.0.146). Tumour-normal sample concordance as well as cross-individual contamination was assessed using Conpair v0.251. Samples were excluded on the basis of quality issues, such as <80% of the targeted regions covered with 20X coverage or more, excessively small library insert size, or cross-individual contamination >5% in either the tumour or normal.

Transcriptome. The FFPE-extracted RNA samples were transcribed, and sequencing libraries prepared using the NEBNext Ultra II Directional RNA Library Prep kit (New England Biolabs) according to the manufacturer’s instructions. Unique dual index tags were applied, and the samples were amplified by PCR using the KAPA HiFi HotStart ReadyMix PCR Kit (Roche) for a minimum of 16 cycles. The libraries were quantified using Accuclear dsDNA Quantitation Kit (Biotium), pooled (8-plex) in an equimolar fashion and hybridised overnight with the SureSelect Human All Exon V5 baits (Agilent). The multiplexed samples were paired-end sequenced using the NovaSeq platform (Illumina) to generate 101bp reads. Reads were aligned using STAR v2.5.0c1552 against the GRCh38 human reference genome alongside Ensembl v103 gene annotation. Expression was assessed by counting reads using HTseq v0.7.253 with the appropriate stranded parameter and subsequently transformed into transcripts per million (TPM). Data quality was assessed by running RNA-SeqQC254 and looking at the total number of counts obtained per sample. We discarded any samples that reported: low expression profiling efficiency <40%, 3 prime bias <0.3 or 3 prime bias >0.5, proportion of reads intersecting rRNAs >2.5, read pairs counted <20×107, had higher number of low quality, ambiguous, alignment not unique than reads counted; or number of genes with five counts or more <14×10355.

Analyses
Whole exome
Somatic variant calling. Somatic point mutations were identified using cgpCaVEMan (v1.15.1)56. SmartPhase (v1.2.1)57 was used to identify MNVs from SNVs called with cgpCaVEMan, while adjacent SNVs were first extracted from VCFs generated by CaVEMan. In addition, indels on the autosomes and chromosomes X and Y were identified using cgpPindel (v.3.10.0)58 (https://github.com/cancerit/cgpPindel). For both cgpCaVEMan and cgpPindel, tumours without a matched normal tissue sample were analysed using a BAM file containing ∼50x exome coverage of simulated 100-bp paired-end sequencing reads with an average 300bp insert size generated from the GRCh38 reference genome sequence, without mismatches.

Gene models and the Variant Effect Predictor (VEP)39 from Ensembl v103 were used to predict the consequences of base changes and indels on proteins. The canonical transcript, as defined in Ensembl, was used to determine the variant consequence. VEP was also used to add custom annotations from the COSMIC (v97)59, ClinVar (update 20230121)40, gnomAD (v3.1.2)60 and dbSNP (v155)61 databases. Common SNVs, defined as variants present in 1% or more of the total population in the gnomAD database or the 1000 Genomes Phase 3 dataset (as indicated by dbSNP v155), were excluded as germline variants. Pindel calls were further refined by retaining only variants with VAF≥0.1 and multinucleotide variants of length 3bp or less. Variants were retained where both the reference allele and alternative allele are ≤25bp in length and both alleles are not >10bp, otherwise the variants were retained only if the tumour VAF >0.25 and sequencing coverage at the variant site was ≥20 in both the tumour and matched normal. Finally, variants within 100bp of exome targeted regions and not flagged by cgpCaVEMan or cgpPindel filters were taken forward for further analysis. For the generation of the oncoplots, the FLAG genes were not included as they are frequently found to bear rare, likely functional variants within the general population62.

Somatic copy number alteration. Somatic copy number alterations (SCNAs) were identified in tumours with matched normal tissue using ASCAT (v3.1.2)63 in all autosomes and the chromosome X. The required hg38 reference files for processing WES data (loci, allele, GC correction and replication timing correction files) were downloaded from https://github.com/VanLoo-lab/ascat/tree/master/ReferenceFiles/WES (git commit ID 29f2fad). The loci and allele files were used as input to the ascat.perpareHTS function, along with tumour and matched normal BAM files. Other parameters used were as follows: genomeVersion=hg38, minCounts=10, min_base_qual=20, min_map_qual=35 and seed=485028101 for reproducibility. A BED file containing the genomic coordinates of the exome pull-down regions sequenced was also provided, as recommended. The sex parameter was determined by the clinical data provided for each patient, and alleleCount (v.4.3.0; https://github.com/cancerit/alleleCount) was also used. The outputs of ascat.prepareHTS were used to run the ascat.correctLogR function with the reference GC and replication timing files, and the ascat.aspcf function was then run with penalty = 70 and seed = 483024451 for reproducibility. Finally, the ascat.runAscat function was run using gamma=1 to estimate purity, ploidy and allele-specific copy number at each loci. Cases with a solution with a goodness-of-fit score <0.9 were considered noisy and excluded from further analysis.

To find significant recurrent amplifications and deletions, the outputs from ASCAT were used to generate input files for GISTIC2 (v2.0.23)64, which requires segment coordinates, the number of markers and a log2-scaled copy number. For each segment obtained from segmentation from ASCAT analysis, the normalised depth log ratio was extracted, and the number of assessed loci within each segment was used as the number of markers. GISTIC2 was run with the following parameters: -refgene hg38.UCSC.add_miR.160920.refgene.mat -gene-gistic 1 -smallmem 1 -broad 1 -brlen 0.75 -conf 0.95 -armpeel 1 -save-gene 1 -gcm extreme -v 20 -ta 0.25 -td 0.25. The hg38.UCSC.add_miR.160920.refgene.mat file is included with GISTIC2. The threshold for the residual q-value for both broad and focal amplifications and deletions was set at 0.10. To remove artefactual recurrent SCNAs, we further removed GISTIC2 wide peaks with an overlap of more than 40% with regions known to be problematic for sequencing and/or read mapping, as defined by the Genome in a Bottle Consortium benchmark union set of all difficult regions (v3.3; https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/genome-stratifications/v3.3/GRCh38@all/Union/GRCh38_alldifficultregions.bed.gz). Finally, we also excluded a wide peak if the list of samples predicted to have a particular SCNA in GISTIC2 did not concur with the original ASCAT results.

Driver genes. To identify driver genes, we used dNdscv v0.1.0 (git commit ID 64f8443)65, which detects genes under positive selection in cancer. The reference gene annotations and covariates files used to run dNdScv were provided in the dNdScv package (RefCDS_human_GRCh38_GencodeV18_recommended.rda and covariates_hg19_hg38_epigenome_pcawg.rda, respectively). The outputs include p-values that were adjusted using the Benjamini & Hochberg method, and the q-value threshold was set at 0.10. A gene was selected if any of the q-values for substitutions and indels combined, substitutions only, missense mutations, nonsense mutations or indels only was below the q-value threshold. For patients with more than one tumour, variants were first merged to create a composite sample by patient prior to running dNdScv.

Hypothesising that the identification of several drivers might have been missed when testing all genes, restricted hypothesis testing (RHT) was performed considering 312 previously reported cancer driver genes17.

OnodriveFML (v2.4.0)66 was also used to identify candidate driver genes. OncodriveFML determines whether the average functional impact score of somatic mutations on an element (in this case, a coding gene) is significantly higher than expected, OncodriveFML was run with the following parameters: genome build, hg38; sequencing, wes; signature method, complement; classifier=CANCER_TYPE; statistic method, amean; sampling, 100,000; sampling max, 1,000,000; sampling chunk, 100; sampling minobs, 10; include indels, true; indel method, max; indel max consecutive, 7. Functional impact scores for the genome were downloaded from Combined Annotation Dependent Depletion (CADD; v1.7) website (https://cadd.gs.washington.edu; https://krishna.gs.washington.edu/download/CADD/v1.7/GRCh38/whole_genome_SNVs.tsv.gz). A regions file consisting of 100bp padded bait set regions was used. A gene was considered significantly mutated if it was mutated in at least 2 samples and had q<0.1.

Mutually exclusive and co-occurring interactions in driver genes (DISCOVER). DISCOVER20 was used to identify mutually exclusive and frequently co-occurring mutations. A binarised alteration matrix was generated for 17657 genes across the 197 primary STs with associated matched normal, where 0 and 1 indicated the absence or presence of a nonsynonymous mutation, respectively. Using the DISCOVER R package v.0.9.4, pairwise tests were performed for all genes which were altered in at least 10 STs. The tests were performed (i) with no stratification by ST subtype, (ii) with stratification by ST subtype, (iii) with stratification by MMR status and (iv) with stratification by ST subtype and MMR status. A false discovery rate (FDR) of 5% was applied to determine statistically significant mutual exclusivities and co-occurrences, focusing on the 28 driver genes called by dNdScv (Supplementary Data 6).

Mutational signatures. To identify somatic mutational signatures for single base substitutions, doublet base substitutions, and indels, somatic mutations from tumours with matched normal tissue were analysed using SigProfilerExtractor (v1.1.21)67. SigProfilerExtractor performs de novo extraction of mutational signatures and assigns known signatures to samples by refitting known COSMIC signatures to the extracted signatures using SigProfilerAssignment68. SigProfilerExtractor was run in exome mode using GRCh38 as the reference and opportunity genome, 500 replicates for non-negative matrix factorisation (NMF) with 1 to 10 signatures. The solution implemented for downstream analysis was the optimal solution provided by SigProfilerExtractor (referred to as the suggested solution). Results from signature decomposition and assignment were taken forward if the cosine similarity between the de novo extracted signatures and signatures reconstructed from assigned COSMIC signatures was ≥0.9. Similarly, on a per-sample basis, signature analysis was considered reliable if the cosine similarity between the original mutational spectrum and reconstructed spectrum was ≥0.9 and there were at least 100 mutations per sample, for each signature type (SBS, DBS, ID).

Germline variant calling. Reads aligned using BWA-MEM50 against the GRCh38 reference, and PCR duplicates were marked using Picard’s MarkDuplicates v2.25.4 (https://broadinstitute.github.io/picard/). Germline variants were identified from whole-exome data from normal samples using GATK v4.2.6.169 following GATK Best Practices70 for cohort calling of SNP and short INDELs. Variants were first called using GATK’s HaplotypeCaller in -ERC GVCF mode with parameters -G StandardAnnotation, -G StandardHCAnnotation, -G AS_StandardAnnotation. Followed by the creation of a genomicsdb database using GenomicsDBImport. Finally, joint genotyping was performed using GenotypeGVCFs. We separated SNP and INDELs using Picard’s SelectVariants v2.27.1. Subsequently, we used Picard’s VariantFiltration v2.27.1 to hard filter variants. For SNPs we filtered out variants: outside 100bp upstream or downstream from the baitset regions, QUAL< 30.0, SOR>3.0, FS>60, MQ<40, MQRankSum<−12.5 and ReadPosRankSum<−8.0. For indels, we filtered out variants: outside 100bp up and downstream from the baitset coordinates used, QD< 2.0, QUAL<30.0, FS>200.0 and ReadPosRankSum <−20.0. High-quality variants functional consequences were annotated using Ensembl Variant Effect Predictor (VEP) v103 using the same parameters, and additional annotations from ClinVar (Updated 2023.11), COSMIC, and gnomAD were added, as the somatic variants. Only indels ≤10bp were kept. To identify variants present within known germline cancer predisposition genes, we looked for variants present within the set of genes reported to be used by England’s National Health Service (NHS) Cancer National Genomic Test Directory v7.2 June 2023 [https://www.england.nhs.uk/wp-content/uploads/2018/08/Cancer-national-genomic-test-directory-version-7.2-June-2023.xlsx] to diagnose cancer predisposition in patients. Only variants with predicted moderate or High Impact over the encoded protein were reported. Variants found in POT1 were removed, as after manual inspection, they appeared to be artefactual.

Ancestry analysis. We projected the germline single-nucleotide polymorphisms from the whole-exome sequencing normal samples into the principal component space of the 1000 Genomes (1000G) project71 using the package AKT v0.3.372. Reported country of birth was aligned with projected super population by visual inspection of 1000G individuals alongside study individuals in the first four principal components, this confirmed that most individuals were either of European or East Asian ancestry.

MMR status. MMR status was defined based on 3 criteria, such as MMR expression (MLH1, MSH2, MSH6 and PMS2) by immunohistochemistry, indels rate and dMMR mutational signatures. Tumours were called dMMR when presenting loss of expression of at least one MMR, high indels rate (>10/Mb) and/or at least one dMMR signature (SBS and/or ID). In contrast, tumours were named pMMR when showing retained expression of MMR, low indels rate (<10/Mb) and/or no dMMR signatures (SBS and/or ID). Of note, tumours with no matched normal presented a higher mutational burden/indels rate due to the lack of accurate SNP filtering.

However, when statistically testing the best predictor for dMMR detection (i.e., indel rate versus tumour mutational burden; TMB (including SNVs, MNVs and indels), only confirmed dMMR and pMMR tumours were selected, in which IHC and somatic/germline variant calling where available. dMMR STs were retained when showing a MMR loss of expression, in addition to a biallelic mutation of an MMR genes.

Transcriptome
Fusion genes. Gene fusion identification was performed using STAR-Fusion v1.10.173, with STAR v2.78a and the Trinity Cancer Transcriptome Analysis Toolkit (CTAT) genome library StarFv1.10 for GRCh38 using GENCODE v37 (Ensembl v103) gene annotation. To assess the coding effect of the fusions and to annotate their presence in previous databases and validate them in silico, we used FusionInspector v2.6.0 by running STAR-Fusion with the parameters —FusionInspector validate —examine_coding_effect —denovo_reconstructFusion. Subsequently, we filtered out gene fusions: fusions with the total number of Junction read counts plus Spanning frag counts <5, FFPM<0.1, not a predicted coding effect, annotated as previously reported in normal tissues of the Genotype-Tissue Expression (GTEx) project or neighbours. Collation of results, filtering, summaries per cohort and plotting was performed using custom scripts. Fusion gene structures were plotted using the Chimeraviz v1.24.074 package in R. In cases where data for multiple samples for the same patient was available, we used only a single sample to report the observed frequency of fused gene pairs within the cohorts.

Cell-fraction estimation of the tumour microenvironment (TME). The Immune Estimation module of TIMER2.0 (Tumour Immune Estimation Resource version 2)75 for immune infiltration estimation was run on the tumour samples. TPM-normalised, gene by sample, expression matrices were provided as input to the TIMER2.0 shiny web server (http://timer.cistrome.org/), selecting AUTO as the cancer type. Using EPIC output from TIMER2, cell-fraction estimation of nine cell types relevant to the TME of skin malignancies were obtained: B cells, cancer associated fibroblasts (CAFs), CD4+ T cells, CD8+ T cells, endothelial cells, tumour cells, and macrophages. These data were used alongside other data for the integrative clustering (described below).

Whole exome and transcriptome
Integrative clustering. Pre-processing of mutation data: SNVs and indels were collated per gene into a binary matrix, with 1s if the gene was mutated in a sample and 0 if the gene was wild-type. The selection of features was based on three criteria: first, genes which were at least mutated in 50 samples were selected. Second, those genes with mutations in at least 15 samples were tested (using a t-test) for association between the presence of mutation and total mutational burden in the tumour. Those genes with a p-value larger than 0.01 after FDR (false discovery rate) correction were selected, as those corresponded to genes where the presence of a mutation is not associated with high mutational burden in the tumour. Finally, based on biological knowledge, cancer genes such as TP53, MSH2, MSH6, MLH1, PMS2 were also included (if not selected by previous rules), and others like TTN, MUC6 and OBSCN were excluded (these genes are FLAGS76, i.e., these are genes often non-pathogenic and passengers, but are frequently mutated in most of the public exome studies). In total, 84 genes were selected.

Pre-processing of RNA-SEQ data: First, transcripts were mapped to gene symbols. In cases where more than one transcript corresponded to the same gene, the transcript with the higher interquartile range (IQR) was selected. Then, the RNA-SEQ read count matrix was converted to transcripts per million kilobases (TPMs) and log transformed (log(1+TPM)). In addition, the voom transformation77 was applied to attain a more suitable distribution for linear modelling. Next, the IQR of each gene was computed, and the 1000 genes that showed the most variability in expression were selected for clustering. We also added to this list the genes selected from the mutation list described above to account for cis-effects, ending up with 1082 genes. This number provided a good balance between high variance in the expression and a reasonable number of features for the computational performance of the Integrative Clustering algorithm. Finally, each gene was standardised to a z-score.

Sample selection: Only one sample replicate per tumour was chosen. Tumours with a diagnosis of SA, SM, SC-E and SC-O were selected, for a total of 170 tumours.

Integrative cluster algorithm. MOFA78 was used to fit a regularised latent variable model-based clustering. This model assumes a series of latent variables to represent k different molecular drivers that predict the observed combination of genomic and transcriptomic values in the tumours. These variables represent continuous activation levels of the different driver processes. The mutation data, represented with a binomial distribution and the expression data, represented with a Gaussian distribution, are connected to the latent driver processes via a parametric joint model. Regularisation is incorporated in the model to reduce model complexity and perform automatic feature selection.

We used the default number of factors (15), studied the robustness of the models fitting 10 runs and comparing the outputs, and selected the best model using the Evidence Lower Bound (ELBO) criteria. We observed that only the first 10 factors explained at least 1% of the variability, and we decided to cluster the tumours based on those first 10 latent variables.

We performed hierarchical clustering using the Ward method on the Euclidean distance matrix and selected the best partition using the silhouette index criteria, composed of 9 clusters. We compared alternative solutions and found that this solution provided also the most meaningful structure.

To interpret the most common features of each group, we compared the frequency of mutations in each cluster and conducted differential expression between the tumours in that cluster and the rest using limma79. Finally, we conducted gene ontology enrichment analysis of those differentially expressed genes using the package clusterProfiler80, correcting the p-values with FDR and selecting those pathways enriched with an FDR <0.05, either over or under represented.

Neoantigen prediction. Class-I and -II haplotyping and neoantigen prediction were performed using the nextNEOpi nextflow pipeline81. Both DNA and RNA were used for predictions where available. Class-I and -II haplotypes for each sample were ascertained using Optitype82 and HLA-HD83, respectively. Variant calling was performed using Mutect284, Varscan285, Manta86 and Strelka287. Variants called by Mutect2 and at least one other caller were selected as candidate neoepitopes. Clonality of neoantigens were determined by ASCAT87 and Sequenza88. Further assessment of peptide-HLA binding was performed by pVACseq from the pVACtools suite89,90. Finally, neoantigenicity was evaluated by a Cauchy-Schwarz index of Neoantigens (CSiN)91.

Pathogen sequence analysis. Screening for pathogens used Kraken292 v2.1.2. Unfiltered paired-end whole exome and transcriptome sequencing FASTQ files were supplied as input, which was run with a 16GB-capped reference database containing RefSeq sequences from archaea, bacteria, viruses, plasmids, humans, UniVec_Core, protozoa, fungi, and plants [https://benlangmead.github.io/aws-indexes/k2] (version of March 2023). The following Kraken2 settings were used to generate i) report-style outputs: —paired —gzip-compressed —use-names —confidence 0.1 —db path/to/DB —report path/to/report —report-minimiser-data —output /dev/null path/to/fastq1 path/to/fastq2; ii) MPA-style outputs: kraken2 —paired —gzip-compressed —use-names —confidence 0.1 —db path/to/DB —report path/to/report —use-mpa-style —output /dev/null path/to/fastq1 path/to/fastq2. A confidence score threshold (—confidence) of 0.1 was applied to compensate the risk of false positives that is inherent to the method [https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown]. Using Kraken2 results, the proportion of minimisers was calculated from Kraken2’s reference database found in each sample for each taxon according to:

Proportion = (Number of distinct minimisers in input read data) / (Total clade level minimisers) (1)

To assess the significance of the results, a normal distribution was used to calculate the right-tailed probability of finding at least m minimisers, out of the total M minimisers available for a taxon, in n reads of a sample, considering the total N reads evaluated by Kraken2 in that sample. The Benjamini-Hochberg method93 was performed to correct for multiple comparisons. Results whose adjusted p-value was smaller than 0.05 were considered statistically significant.

Impact of alterations in pathway activations. Hippo and Yap1 pathway activations were computed based on the REACTOME definitions downloaded from the Molecular Signatures Database94 using Gene Set Enrichment Analysis with the R package gsva95. Next, we identified potential molecular aberrations that might have an effect on the inactivation of these pathways96. Each mutation was classified as follows: for oncogenes, missense variants, splice acceptor variants and splice donor variants were classified as activating mutations. For tumour suppressor genes, frameshift variants, start-lost, stop-gained and stop-lost variants were classified as truncating mutations. We also computed loss of heterozygosity using the copy number data. In order to test the association of the presence of these mutations with changes in the activation levels of the pathway, we performed a linear model and conducted a t-test on the parameters of the model.

MLH1 promoter hypermethylation testing. Analysis performed using tumour FFPE tissue derived DNA, which was bisulphite converted using the EpiTect Fast Bisulfit Kit. The methylation status of the tumour tissue was assessed by methylation sensitive high resolution melting analysis (MS-HRM), using MethylDetect MLH1 primer kits and EpiTech HRM PCR master mix on a Rotor-Gene Q 5Plex HRM platform. The MethylDetect MLH1 kit contains primers that amplify a region covering 5 CpG sites in the 5’ UTR region of the MLH1 gene (NM_000249.3), along with high-, low- and non-methylated control materials. MS-HRM data analysis was conducted using the Rotor-Gene Q Software v2.3.4 (build 3).

Ethics approval
Ethical approval was obtained from the different Institutes contributing to this project. Genetic analysis was further approved by the Sanger Human Materials and Data Management Committee.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is a large and comprehensive study of genomic alterations in sebaceous neoplasms, reaching from benign to malignant tumors. The data generated represents a considerable resource, which is important for the scientific community. While the scope of the study and the number of samples is impressive, there are several significant weaknesses. Overall, it is not fully clear what novel aspects of these tumors has been revealed. The 5 mutational patterns are not fully convincing and the analysis of copy number changes and structural rearrangements not compelling as presented.\nThe manuscript is difficult to read with many nested sentences full of abbreviations and not always clear meaning and partially problematic grammar. More importantly, no coherent picture is distilled by integrating the multi-omic data to extend the knowledge on the biology of these tumors, i.e. what signaling pathways are consistently altered, what constitutes the barrier between benign and malignant, and how any of these features might relate to clinical course. Outcome data is mentioned to exist but does not seem to feature in the analysis.\nSpecific points:\n1) Figure 1B y-axis labels should be clarified. What additional tumors qualify? Other STs. Could be phrased more clearly and succinctly.\n2) Figure 1C is not referenced in the text and the labels are not clear. The pictures are clinical and histological and not clinicopathological\n3) It would be good to show a lollipop diagram for RREB1 as well. Is there any information that can be gleaned about the functional relationship between RREB1 and HRAS mutations? For example, are particular types of RREB1 mutations associated with HRAS mutations? What are the anticipated functional consequences of RREB1 mutations?\n4) The last sentence on mutually exclusivity is difficult to parse and the pattern of exclusivity not apparent (to me) from figure 2. It looks like many tumors with TP53 mutations also have mutations in RREB1, HRAS, or NOTCH1\n5) The data on copy number changes implicating individual genes is not compelling. GISTIC scores in fig 5 looks noisy rather than showing bona fide copy number changes. This is particularly true for SA and SM, which seem to have overall ‘flat’ genomes. The sentence on page 8 line 342 onwards on amplifications of specific genes cannot be understood.\nAlso, what constitutes a ‘significantly focal amplification’? Compelling examples of the raw copy number data for example for HRNR and KIT amplifications should be shown for this section to be more compelling, ideally with additional supporting information nominating these genes (FISH, IHC). A better integration with mutation data is desirable when nominating specific genes based on copy number change. E.g. how would an amplification of GNAS be plausible in this tumor type, in which no activating mutations of GNAS are seemingly present?\nThe wording and grammar in these paragraphs is awkward e.g. ‘In contrast, SC-Es showed significantly chromosomes 1 and 20 amplifications, and significantly 17p deletions, which includes TP53 (7/36 tumours, of which 6 also had a somatic TP53 mutations’\n6) Additional consideration and discussion should be given to the PAK2 fusions. The fusions with PAK2 at the 5’ end appear to be novel/unprecedented. Is there a reciprocal fusion product that was perhaps missed?\n7) The finding of an identical fusion of RCOR1::GRHL2 needs additional discussion as it implies a clonal relationship between the two SAs of this patient.\n8) The finding on CIN primarily in pMMR tumors is interesting and mirrors the findings in CRC. This indicates, as you state, that genomic instability is necessary for the malignant transformation. Any insights on why and how?", "reviewer_2": "General comments:\nIn this multi-center international study, Ferreira et al. present the largest exome and transcriptome analysis of sebaceous tumors (ST) to date. This allowed investigators to identify mutated genes under positive selection (drivers) for the first time in this tumor type, which provides an incredibly valuable resource to the cancer research community. In my opinion, this is an extremely important study for the ST and cancer genomics fields.\nThe initial data processing, including variant calling, mutation signature extraction, and copy number profiling, was conducted using well-established software and pipelines. However, several downstream analyses are problematic, and the presentation of the results can be significantly improved to enhance clarity and better guide the reader. Furthermore, in my opinion, parts of the neoantigen and ligand-receptor analyses are too speculative, and they lack additional data or experimentation for validation. These sections should be modified for clarity and rigor, and their results presented more cautiously or the authors may want to consider omitting these analyses from the study.\nMy recommendation is to improve the whole exome and transcriptome analysis as outlined in detailed below and re- consider inclusion of the ligand-receptor analysis and parts of the neoantigen analysis.\nConcerns/suggested revisions:\nSomatic landscape of sebaceous tumours\nDriver Identification Analysis:\n[1] The paper reports significantly mutated genes based on a single software tool (dNdScv), without the necessary validation steps (e.g. filtering for genes that are not expressed or have low VAF), raising concerns about false positives that are ubiquitous in the context of high TMB tumors. (For example, SPANXB1 is reported as a driver in Figure 2, but most mutations in this gene appear to have a suspiciously low VAF). To ensure robustness, the authors should filter their results to exclude likely false positives and, importantly, provide additional supporting evidence for the novel candidate drivers:\n● What is the variant allele frequency (VAF) distribution for each gene?\n● Do these variants mostly originate from a specific signature or positional site?\n● What proportion of variants are seen in the RNA-Seq data for each gene?\n● Are there significant co-occurrences of loss-of-function mutations with loss of heterozygosity (LOH) for the putative tumor suppressor genes (TSGs)?\n[1.1] The authors of the dNdScv software have excluded hypermutator samples (>500 mutations) from their main analysis (iMartincorena et al, 2017, PMID:29056346). Such precautions are not mentioned in the method section of the current manuscript. Did they authors try running dNdScv exclusively on tumors with fewer than, say, 1000 or 2000 mutations to see how it affects the results?\n[1.2] There are 24 samples with no matched normal. It is important to exclude these samples from the dNdScv analysis as they could contain germline variants misclassified as somatic (I’ve noticed VAF modes near 0.5 for many of these samples). Could the authors specify which samples have been included in their analyses, along with any exclusion criteria?\n[1.3] Possible additional figures for the driver analysis (in main or supplementary): A barplot showing the relative proportion of each variant class (missense, frameshift stop gain etc) in each gene or selected genes. A scatter plot showing dN/dS:missense vs dN/dS:Lof similar to Figure 2B in Martincorena et al.(2017, PMID: 29056346). These figures would help support the analysis results and importantly inform about the action mode of the candidate drivers (i.e. oncogene vs.TSG).\n[1.4] It would be helpful if the authors cross-reference their gene list with the results of some pan-cancer analyses, such as Bailey et al. (2018, PMID: 29625053) or Martincorena et al. (2017, PMID: 29056346). This information could be easily summarized in a supplementary figure having genes as rows and cancer types as columns, with each tile colored according to their report status (e.g. oncogene, potential oncogene, tsg and potential tsg, which is available in Supplementary Table 1 of Bailey et al.). It would be useful for the readers for genes that haven’t been reported as significantly mutated, what does the literature say about them that could support their role as an oncogene or tumor suppressor?\n[1.5] I speculate the authors may have missed the identification of several drivers by testing all genes as they did impose a burden on significance levels after multiple testing corrections. In addition to their global analysis, I think it may be worth performing a restricted hypothesis testing (RHT), such as in Lawrence et al. (2014, PMID: 24390350) by considering only previously reported genes, for example those from Bailey et al. (2018), and those reported as drivers in sebaceous tumors. This additional analysis would be straightforward, requiring only the subsetting of the global results to these genes before applying multiple testing corrections. This may be valuable to the readers.\n[1.6] The main oncoplot could be enhanced in several ways.\n• Adding copy number track for the most frequently affected genes (deleted TSG or amplified oncogenes).\n• Some drivers with less supporting evidence or of lower interest could be removed from the main figure for simplicity and clarity (but shown in a full supplementary oncoplot).\n• Annotation tracks with tumor purity, RNA subtypes and mutational signature groups could be added.\n• Mode of action of each driver (TSG vs oncogene) could be highlighted.\n[1.7] The authors only report the gene names with their corrected p-values in the supplementary material. It is important to include the complete dNdScv output in the supplementary material to provide a comprehensive view of the analysis and facilitate further validation and reproducibility of results.\n[1.8] On lines 199-200, the authors write, \"Interestingly, when considering the ST subtypes individually, HRAS was only identified as a driver gene in the benign tumors.\" What is the uncorrected p-value for HRAS in the SC-E and SC-O subgroup analyses? Would it pass RHT? Even if not significant, the oncoplot shows several HRAS mutations in SC-E and SC-O, most of which are known activating mutations at G12/13 and Q61. Moreover, HRAS harbors no silent mutations in these tumors, strongly suggesting it is an actual driver, albeit less frequently mutated. Accordingly, it would be more appropriate to drop this statement or at least revise it to \"HRAS was not reported as significant in (...) probably due to its lower mutation frequency and the smaller size of these groups (...),\" and mention the presence of these hotspot mutations. In general, the absence of a significant result for a particular gene in a subgroup analysis cannot be used on its own to suggest an absence of “driver” activity.\n[1.9] In the discussion, lines 491-492, the authors state: “Use of somatic mutational profiling may help in the differential diagnosis of benign versus malignant STs, which can be challenging.” It would be helpful if more clinical context could be provided here, particularly around the standard histopathological diagnostic workflow. What are the gold-standard IHC stains, histopathological correlates, or clinical decision algorithms currently used to differentiate between the four types of sebaceous tumors? A deeper discussion on how mutational testing could provide clinical input would be ideal, given the clinical correlation angle in the discussion. Figure 2 does not clearly identify a mutational profile with high sensitivity or specificity for diagnosis. Could the authors consider using classifiers (e.g., Random Forests) and measure their performance with cross-validation?\nMutually exclusive and co-occurring mutations:\n[2] The “mutual exclusions and co-occurring interactions” analysis is problematic as it combines the results of a Fisher’s test with the output of DISCOVER, a tool specifically developed to replace the inappropriate usage of the Fisher’s test in such analyses. The Fisher’s test does not control for TMB, leading to a high number of false positive co-occurring interactions (as shown in Figure 2B). In contrast, DISCOVER accounts for TMB in each sample and gene, thus limiting these spurious results. Combining the results from both approaches does not make sense, as exemplified by the contradictions between the interactions reported in Figure 2B and the DISCOVER analysis output provided in Supplementary Table 6 (e.g., significant co-occurrence between ZNF750 and BICRA while Supplementary Table 6 reports significant mutual exclusion). I would suggest that the authors should rely solely on DISCOVER (or similar approaches controlling for TMB), and completely exclude the Fisher’s test from this analysis.\nMutational signatures of sebaceous tumours:\n[3] The mutational signature analysis was not comprehensively explored and presented, and a global view of the cohort is lacking, specifically related to the relative contribution of the main signatures in each sample. The presentation of the results in Figure S5 show too many signatures simultaneously, and additional figures and analyses would be necessary to support the claims and observations made in the main text and discussion. In my opinion, an improved global view of the mutational signatures would warrant inclusion as a main figure.\n[3.1] Figure S5 needs clarification. The graphical legend mentions */PD A, */PD NA, and */PD NV in three different colors, but none of these are explained in the textual legend. There is a reference to the blue and red labeling of the sample names, but this is confusing since \"blue\" samples (no significant signature) actually show some signatures. Additionally, the legend states that \"one in red text was not included,\" but half of the samples are in red.\n[3.2] The method section states that SigProfilerExtractor was used for de novo extraction, considering solutions with 1 to 10 signatures by NMF. However, Supplementary Figure 5 shows 21 signatures (17 SBS, 3 ID, and 1 DBS), and Supplementary Table 8 reports mutation counts for 27 signatures. This discrepancy needs clarification. It is essential to explain the criteria for determining the final set of signatures and how they relate to the initial NMF analysis to ensure consistency and transparency.\n[3.3] I strongly suggest the authors include a main figure summarizing the mutational signature analysis. For example, a clustered heatmap of each signature log10 counts (as rows) in each sample (columns) clearly shows the 3 main mutational profiles within this cohort: MMR (SBS15 + ID2 + ID7 + high SBS1 + high SBS5), UV (DBS1 + SBS7a + SBS7b + low SBS1 + low SBS5) and others (low SBS1 + low SBS5). Annotation tracks for samples could include the MMR status, log10TMB, diagnosis and presence/absence of LoF mutations in MSH2 & MSH6. An annotation track for the signatures could inform the reader about their aetiology (UV, MMR etc). Such heatmap can be easily produced with R’s pheatmap package and the ward.D2 clustering method. To make it more concise and suitable for a main figure, some highly correlated or same aetiology signatures could be merged (e.g. SBS7a & SBS7b) and low contributing signatures could be dropped and deferred to supplementary figures (if necessary). Alternatively, the authors could use the W matrix from NMF or apply PCA on the 96-SBS mutation proportion matrix to create a scatter plot of all samples colored by cancer types, MMR status and other variables of interest.\n[3.4] Lines 268-269, the authors write, “As expected, in the dMMR cases, the high TMB was attributed to dMMR-associated mutational signatures.” What proportion of TMB is attributed to dMMR signatures? Could the authors please include a figure supporting this? It seems that other signatures (SBS1 & SBS5) contribute more to TMB than dMMR-associated signatures, even in dMMR samples.\n[3.5] Lines 163-166, the authors write, “Interestingly, indel rate rather than TMB, more accurately predicted MMR status (indel rate p-value < 2.2e-16 versus TMB rate p-value = 2.203e-08, Welch two samples t-test)”. I have a few concerns regarding this analysis. First, the considered data doesn’t meet this test's assumptions about normality. Second, p-values are not an appropriate measure of prediction accuracy. Third, according to the methods section, MMR status was in part determined by the indel rate (lines 772-774), which means using the indel rate to predict MMR status is circular reasoning. These issues should be addressed, for example by comparing ROC-AUCs and and an MMR classification that doesn’t rely on indel rate or TMB.\n[3.6] On lines 501-513, the authors write “this study highlights the importance of recording the precise location of SC-O tumours, and specifically restricting use of the term SC-O for only those arising on the eyelid. Indeed, the sample coming from the canthus showed a different phenotype from the other SC-O, being characterised by UV mutational signatures rather than being pauci-mutated”. Could the authors please include a figure and additional analyses or references supporting this claim, for example by showing differences in mutated drivers or gene expression that suggest it is more similar to SC-E than SC-O?\n[3.7] Lines 1377-1379: “SBS” should be “DBS”. Also, a typo in “mutationa”.\n[3.8] In some other cancers, MMR status is linked to survival (e.g. Endometrial). Did the authors investigate whether in this cohort, MMR status is linked to survival in SC regardless of subtype? Authors note that sample size in SC was too small to measure but it may be helpful for readers to see this information anyways even if not statistically significant.\nIntegrative cluster analysis:\n[4] The results of this analysis are rather cryptic, and the methodological details raise concerns about their validity. First, the expression data was combined with the mutation profiles of approximately 80 genes selected on the basis that they are frequently mutated in this cohort. Given that a typical tumor contains fewer than 10-15 driver mutations (referring here to Martincorena et al, 2017) and that these sebaceous tumors have huge TMBs (nearly half of them have a thousand mutations or more), we can be sure that most of the mutations in these 80 genes are passenger mutations. While a few relevant genes such as TP53 and MSH2 have been manually added and others like TTN have been manually removed, Supplementary Figure 4 suggests that the limited number of potentially meaningful mutations is still overwhelmed by the number of passengers.\nSecond, the method contains many steps with somehow arbitrary choices (e.g., a first selection step out of 16 fitted models followed by a second selection step for the number of k-means clusters), making it very difficult to assess the robustness of the findings. What happens if we use another fitted model (many had a similar balance of variance explained and a low number of clusters)? What happens if a different k is chosen or if we change the gene sets?\nI should note, even if the authors are able to address the concerns described above, and although I appreciate these exploratory types of integrative analysis, I would argue they cannot replace the more well-established analyses that many readers would expect, if only for clarity and easier interpretation of the results. In other words:\n[4.1] How many (purely) transcriptional groups seem to be there? If the exact number is hard to determine, the smallest number that make sense (visually) should be the starting point.\n[4.2] Once these groups have been identified, do they correspond to the ST subtypes (SA vs. SC-E vs. SC-O)?\n[4.3] Are they enriched for specific mutations or CNAs?\n[4.4] Which immune and other non -tumor cell types are present in each of these bulk tumors? How does this correlate with tumor purity?\n[4.5] Which pathways or genes seems differentially expressed between these groups? Could this be confounded by tumor purity or non-tumor cells or is there evidence that they could represent genuine differences between the cancerous or neoplastic cells?\nClinical data:\n[5] Supplementary Table 1 is a somewhat ambiguous and includes a mix of patient-related and sample-related data with several issues. I recommend that the authors provide two separate tables: one for the patient-related data and one for the sample-related data, each with appropriate IDs (patient and sample).\n[5.1] Some fields lack sufficient description in the legend, such as \"Age,\" which does not specify whether it refers to the age at diagnosis or at sample procurement. Similarly, “Recurrence”, “TimeAfterInitialDiagnosis_Rec(months)”, “Metastasis” and “TimeAfterInitialDiagnosis_Met(months)” also needs some explanation.\n[5.2] Column AE has no name. It contains “Father, Brother” for sample “PD42477a”, but “FamilyCancerLynch” says “no” for this sample.\n[5.3] Some field names differ between the table and the legend. For example, \"MMRstatus\" is mentioned in the legend, but \"MMR.d.MMR.p_simplified\" appears in the table. Same for some IHC fields.\n[5.4] The MMR status is not written correctly for PD42171a (MMR-p instead of pMMR).\n[5.5] For \"MLH1Hypermathylation,\" a mix of blank and \"not tested\" is used. Please use a unique value (either blank or \"not tested\") for all cases to maintain consistency.\n[5.6] Typos: “Hypermathylation”, “Supplemntary”\n[5.7] There are a number of issues in this table. It would be ideal to ensure that all fields used in the analyses are clearly documented and are error-free. I would recommend fields that cannot be used in any reasonable analyses, because of incomplete or ambiguous data, should be dropped.\nGermline alleles in sebaceous tumour patients:\n[6.1] This section could be further developed. Associated Figure 6 contains 27 genes, but only a handful of these genes are mentioned in the text. I would recommend improving figure (e.g., removing genes not discussed and adding other plots showing how mutations correlate with phenotypes) or moving it to a Supplementary Figure.\n[6.2] VAF and other important information is missing from Supp. Table 12. While the method section states that germline variants were identified with GATK HaplotypeCaller and GenotypeGVCFs, Supp. Table 12 includes fields such as “Tumor_Seq_Allele2” which are typically associated with somatic variant calling. Could the authors clarify this?\nCopy number landscape of sebaceous tumours:\n[7.1] This section could be better developed. It would be ideal if both Figures 4 and 5 could be simplified to include only the key panels that are clearly highlighted to emphasize points made in the main text.\nNeoantigen and ICI response prediction:\n[8] The last section of the paper (lines 413-482) reports results based on a number of software tools (CODEFACS, TIMER2.0, EPIC, LIRICS, CIBERSORT, pVACseq, and others) and several external datasets that I evaluate as highly speculative. There are insufficient details provided to properly understand the methodology and the results, which are supported by a single, not particularly informative, main figure. At this stage, I would recommend removing this section to ensure clarity and rigor of the whole exome and transcriptome analysis of the paper, or significantly revise this section with more clear details on how the analysis was performed. However, I should note I anticipate the results will likely still be speculative and omission of the ligand-receptor analysis should be considered if validation by IHC/IF staining studies are not being included.\nOther:\n[9.1] There is some inconsistent naming and abbreviation usage for the tumor subtypes across the main text, figures, and their legends. For instance, \"SM\" is used for sebaceoma in Figure 1 and Table S1, while only \"S\" is used in Figure 3 and S5. Similarly, \"PeriOcular,\" \"Peri-ocular,\" and \"ocular\" are all used to refer to ocular sebaceous carcinoma. Please standardize the terminology and abbreviations for each disease subtype.\n[9.2] According to ST1, sample PD43886c and PD48829c have been used in DNA analysis (DNA_Analysis == Yes). However, I couldn’t find any mutations associated with these samples in ST4.\n[9.3] SF1A has no legend describing the meaning of each color and does not seem appropriate to clearly support the claims.\n[9.4] The cohort size is relatively large, considering the rarity of these tumors. However the SC sub-group is quite small, particularly the SC-O subgroup which may limit the power of the analyses. More discussion of this (and other limitations) would be appreciated in a short, dedicated section.\n[9.5] Suggestion: Is it necessary to note that no recurrence or metastasis arose in SA or SM (Line.149) as these are benign tumors that do not usually (or at all?) metastasize. If they do, this should be noted in the text.\n[9.6] Suggestion: It would be ideal to list the 6 countries that these tumors come from (line 81) and mention in discussion any potential biases from tumor selection based on ethnicity or geographical location – for example, authors mention that SC-O is more common in people of Asian descent (40-60% of eyelid tumors). Have the authors investigated the ethnic makeup of their cohort? Is it representative of the global or regional prevalence/distribution of these tumor subtypes? If not, how might this affect the results of their study?", "reviewer_3": "Authors investigated the genetics and mutational signatures of sebaceous tumors and reached into the conclusion that there are multiple subtypes present. The study also has an interesting observation that some of these tumors also have lynch signatures. The lynch signatures should be investigated further using current non tumor and historical patient samples to reach a definitive conclusion that these mutations were not acquired later in development. Utilizing some published material to support Sebaceous tumor prevalence in Lynch syndrome would also be informative. Authors utilize available pipelines to perform a lot of their computational analysis, it would be helpful if same pipelines are used on some examples of published historical skin tumors such as Melanoma, squamous carcinoma, UV mediated Merkel cell carcinoma for comparison in these figures such that the reader could make comparisons of mutational burden etc. these results may be skewed by the pipelines used. Authors make some suggestions about the immunogenicity of RCC tumors, this is still currently a puzzle since we don't know what makes RCC tumors highly immunogenic. First the Mutational burden for sure does not explain RCC immunogenicity and secondly the CD8 T cell infiltration in RCC tumors is actually deleterious. So the discussion and comparison to these tumors would need to be carefully worded and considered in the text. Authors make conclusions about the immunogenicity of the tumors. This part requires additional work. First the authors again utilize some pipeline to do their analysis, this would be valid if there are some controls presented for comparison. So the Neoantigen score for Melanomas, squamous carcinomas and UV induced merkel cell carcinoma utilizing published data should be presented as a control comparison. Additionally the authors present histology data in figure 1. As a supplement an immune score calculation utilizing CD4/CD8 infiltration should be presented for each tumor since the authors have access to FFPE blocks. This would make the interpretation of the all neoantigen loads the authors are presenting and would have to validate some of the observations"}, "author_response": "We thank Reviewer1 for their constructive comments and insights. As mentioned, we have collected the largest ever dataset of sebaceous tumours, comparing benign (sebaceous adenoma and sebaceoma) and malignant (extra-ocular and peri-ocular sebaceous carcinoma) cohorts which we have examined through extensive computational analyses.\nWe were able to identify several novel aspects in this study:\n1. We discovered (and validate statistically) new driver genes (for example HRAS, MSH2/MLH1, RB1) not previously reported in sebaceous tumours and profiled the full spectrum of mechanisms by which previously reported genes are altered to become drivers (i.e. somatic mutation and/or somatic copy number alterations (deletion, amplification or cnLOH)). For example, we detected copy neutral LOH of mismatch repair (MMR) genes and amplification of other drivers such as GNAS, PLAG1 and Hornerin. We also describe new fusions and importantly show that PAK fusions are not restricted to apocrine poromas – this is of clinical/diagnostic significance. Remarkably, we also identified pathogenic (exonuclease proofreading defects) POLE/POLD1 somatic mutations, not previously reported to be linked to sebaceous tumours.\n2. We also show, for the first time, that tumours may take different paths to tumourigenesis. For example, Lynch patients can carry STs with a germline and somatic mutations, as well as a germline mutation and a cnLOH event. While other tumours which are pMMR, specifically SCs, exhibited biallelic inactivation of TP53 co-occuring with ZNF750 and/or RB1 disruption (particularly SC-O), suggesting TP53 might accelerate ZNF750/RB1-driven SC as occurs for BRAF p.V600E in melanoma. In support of this hypothesis, concomitant TP53/RB1 mutations were identified in some metastatic SC cases. This study reinforces the important role of RB1 in malignant STs, with inactivation occurring via somatic mutation, focal deletion (in SC-E) and 13q deletion (in SC-O). Moreover, the spectrum of the somatic mutations in the tumours across the same individual can vary (see patient PD43887).\n3. We clearly found that sebaceous tumours have a propensity to develop in the context of a high mutational burden including mismatch repair (MMR) deficiency and UV-damage as previously reported, but also POLE/POLD1 somatic mutation and APOBEC activation. In addition, MLH1 promoter hypermethylation was revealed as a rare event in sebaceous tumours, associated with somatic dMMR.\n4. We also show that clinical variables such as developing a sebaceous tumour outside the head and neck area, having multiple sebaceous tumours, and presenting with another non-sebaceous Lynch syndrome related neoplasm are good Lynch syndrome predictors. In addition, histopathological variables such as no epidermal connection, tumour delimitated by a collarette, KA-like architecture and diagnosis of SM are also good LS predictors, with clinical variables being better predictors.\n5. We demonstrated that indel rate is more reliable than TMB in defining dMMR. In skin, UV-damage increases the TMB and the level of these mutations differs depending on the site and age of the patient, while a high indel burden is largely acquired as a result of dMMR.\n6. The method applied to detect Lynch syndrome via the analysis of colon tumours cannot be directly applied to sebaceous tumours. Indeed, unlike colorectal cancers, somatic BRAF p.V600E and MLH1 promoter hypermethylation testing is not relevant to STs and cannot be employed to distinguish sporadic from germline dMMR tumours owing to the absence of BRAF p.V600E mutations in these tumours and our observation that hypermethylation of the MLH1 promoter is a rare event in sporadic dMMR STs.\n7. The use of somatic mutational profiling may help in the differential diagnosis of benign versus malignant STs, which can be challenging. For example, malignant STs were driven by chromosomal instability and typically showed biallelic TP53 inactivation, as well as TP53 mutations co-occurring with ZNF750 and/or RB1 mutations. In addition, dMMR SC-E was characterised by mutations in NOTCH1, RREB1, and DRD5, whereas SC-O was characterised by significant focal amplifications encompassing GNAS/SALL4 and LYN/PLAG1, significant broad 13q deletions (encompassing BRCA2 and RB1), and 8q amplifications (encompassing MYC), with amplification of the MYC proto-oncogene previously reported by targeted sequencing and FISH11. Lastly, both SC-E and SC-O showed significant amplification of chromosomes 1q21.3 and 20. The former amplification encompasses HNRN (Hornerin), which is part of the epidermal differentiation complex on 1q21. Hornerin is expressed by the ocular surface (corneal and conjunctival aspects, including meibomian and Zeis glands), the lacrimal apparatus and the skin, and has antimicrobial and protective (keratinisation) functions 47,48, reinforcing its potential role in sebaceous carcinoma. STs, with the exception of SC-O, can be the first sentinel sign of LS and 3 clinical criteria and 4 histopathological criteria – including (i) the presence of multiple STs, (ii) STs being located outside the head and neck area, (iii) the presence of a LS-associated cancer, in the clinical model, and (i) the absence of epidermal connection, (ii) the delimitation of the ST by a collarette, (iii) a KA-like architecture, and (iv) the diagnosis of SM in the histopathological model, together with loss of MMR expression by IHC – were found to significantly associate with a LS diagnosis (the clinical model being a better predictor) and thus should raise the index of suspicion that could provoke a referral for germline genetic testing. Of note, while the idea of generating a classifier is one we feel is very validate, a comparable dataset of STs to the one we describe in our paper does not exist which makes validating such a tool problematic and thus we feel this analysis is premature. Further, this study highlights the importance of recording the precise location of SC-O tumours, and specifically restricting use of the term SC-O to those lesions arising on the eyelid. Indeed, the sample coming from the canthus showed a different phenotype from the other SC-O and was characterised by a UV mutational signature, rather than being pauci-mutated, the pattern observed in SC-Os.\nDespite the enormous size of this cohort and the deep characterization we are still not powered enough to make conclusions about patient outcome because at the time of writing only 4 patients had died of their disease.\nAddressing specifically the point about multi-omic/data integration – in Figure 3 we integrate our SNV and expression catalogue to define 9 statistically significant groups.\nRegarding the grammar, we have reviewed this carefully and had a professional editor review the text prior to resubmission.\nSpecific points:\n1) Figure 1B y-axis labels should be clarified. What additional tumors qualify? Other STs. Could be phrased more clearly and succinctly.\nWe thank Reviewer1 for making this interesting comment. The figure shows all additional tumours in each patient that were not sebaceous tumours. The y axis was amended for clarity.\nThe designation of Lynch syndrome related cancers comes from the “Revised guidelines for the clinical management of Lynch syndrome (HNPCC): recommendations by a group of European experts” guidelines/paper published in Gut. 2013 Jun;62(6):812-23 (https://gut.bmj.com/content/62/6/812). We have added the reference into the Figure1 legend (please see revised Manuscript). As discussed in the guidelines there remains some ambiguity as to which tumours are/aren’t Lynch related but we followed these consensus guidelines which are used clinically and thus represent the state of the art at this time.\n2) Figure 1C is not referenced in the text and the labels are not clear. The pictures are clinical and histological and not clinicopathological\nWe thank Reviewer1 for pointing out this mistake and we have amended the text and Figure 1 accordingly (please see revised Manuscript and Figure 1). We also added Figure 1C to the Supplementary Discussion with more detailed labels and an explanation of the histology.\n3) It would be good to show a lollipop diagram for RREB1 as well. Is there any information that can be gleaned about the functional relationship between RREB1 and HRAS mutations? For example, are particular types of RREB1 mutations associated with HRAS mutations? What are the anticipated functional consequences of RREB1 mutations?\nWe thank Reviewer1 for this suggestion and we added a RREB1 lollipop diagram in Supplementary Figure 6.\nFocusing first on the highly recurrent RREB1 frameshift (p.E1502Rfs*81) mutation, we found that it was not significantly co-occurring with HRAS mutation (Chisq-test: p-value=0.624, or Fisher test: p-value=0.602).\nRREB1\nWT frameshift_p.E1502Rfs*81\nWT 92 24\nHRAS\nMut 61 20\nIn the same way we do not find an obvious association with other events, for example when we examine significance using all samples with frameshift events.\nRegarding the function of the RREB1 mutations – RREB1 is involved in regulation of cell proliferation, transcriptional regulation and DNA damage repair. We are also aware of unpublished data suggesting a role for mutation of this gene in immune evasion. From the mutations we called, other than loss of RREB1 being tumour promoting, we are unable to glean any obvious mechanism. It is possible/likely that loss of RREB1 contributes to multiple mechanisms of tumorigenesis, and this is something to be explored in future studies.\n4) The last sentence on mutually exclusivity is difficult to parse and the pattern of exclusivity not apparent (to me) from figure 2. It looks like many tumors with TP53 mutations also have mutations in RREB1, HRAS, or NOTCH1\nTile plots do not consider gene length or sequence context or variant consequence, which is considered in the statistical model deployed within DISCOVER meaning that it can be difficult to decern genetic interactions by looking at these profiles alone. We note that our analysis with DISCOVER (and also with Fishers; removed from the revised manuscript as requested by the reviewers) was performed without stratification for tumour type, thus including the SC-O cases where the mutual exclusivity between the abovementioned genes is more evident.\nWe have reviewed these associations and they remain statistically robust.\nTP53 and HRAS were statistically mutually exclusive, which was demonstrated by DISCOVER (adjusted pvalue = 0.00091) and Fisher (one-tailed pvalue = 0.0381) tests. Although we have samples with both mutations, it was defined as statistically significant, likely driven by SC-O.\nWe also ran DIALECT (Shuaibi A et al. 2024; https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11071465/), an alternative algorithm for identifying co-occurring and mutually exclusive interactions. However, even after corresponding with the authors and rebuilding their code we were unable to generate meaningful output. We therefore conclude that DISCOVER remains the state-of-the art for such analyses and thus we report these results in the paper.\n5) The data on copy number changes implicating individual genes is not compelling. GISTIC scores in fig 5 looks noisy rather than showing bona fide copy number changes. This is particularly true for SA and SM, which seem to have overall ‘flat’ genomes.\nWe have re-run the analysis as shown below.\nFigure R1/Figure 6 – Significant focal somatic copy number alterations in sebaceous tumours. Representation of the GISTIC2.0-identified by ST subtype significant focal peak regions (≤ half a chromosome arm) in (A) sebaceous adenomas (SA), (B) sebaceomas (SM), (C) extra-ocular sebaceous carcinomas (SC-E) and (D) peri-ocular sebaceous carcinomas (SC-O). Chromosomal location of the peaks is indicated by the cytoband, with amplifications (gains) shown in red and deletions (losses) shown in blue. The presence of known tumour suppressor genes or oncogenes within these peak regions are labelled. For the G-score, positive values refer to gains, while negative values are related to deletions. (E) Oncoplot showing the individual tumour samples with the amplification and deletion regions shown in A-D. Samples with no significant gains or losses are indicated by ‘NS’. Five samples were not included in the GISTIC2.0 analysis (see Methods), so no data is shown (indicated by ‘No’). Abbreviations: MMR, mismatch repair; dMMR, deficient in MMR; No, not analysed; NS, not significant; pMMR, proficient in MMR; SA, sebaceous adenoma; SM, sebaceoma; SC-E, extra-ocular sebaceous carcinoma; SC-O, peri-ocular sebaceous carcinoma.\nNote: *focal peaks were considered as significant when reaching the 4 criteria: 1. residual q-value < 0.1, 2. peak size >= 100,000bp, 3. concordance with ASCAT CN calls ≥ 0.75, 4. portion of peak overlap with Genome In A Bottle ‘difficult regions’ <= 0.4 (Supplementary Table 10); **approximate cytobands are shown with related genomic coordinates available in Supplementary Table 10.\nThe sentence on page 8 line 342 onwards on amplifications of specific genes cannot be understood.\nSC-Os showed the greatest number of focal SCNAs (Figure R1/Figure 6). Oncogenes in these regions included GNAS, NFATC2, PTK6 and SALL4 (14/15), and LYN and PLAG1 (14/15). We also observed deletions of the TSGs EZH2, TRIM24, CNTNAP2 and KMT2C (5/15). Of note, the focal SCNAs in SC-Es only included known TSGs, specifically FOXO1 and RB1 (8/36), and ARHGEF10 (9/36).\nAlso, what constitutes a ‘significantly focal amplification’?\nWe have changed significantly focal to significant focal amplification.\nTo find significant recurrent amplifications and deletions, the output of ASCAT was used to generate input files for GISTIC2 (v2.0.23)65, which requires segment coordinates, the number of markers and a log2-scaled copy number for each segment. GISTIC2 was run with the following parameters: -refgene hg38.UCSC.add_miR.160920.refgene.mat -genegistic 1 -smallmem 1 -broad 1 -brlen 0.75 -conf 0.95 -armpeel 1 -savegene 1 -gcm extreme -v 20 -ta 0.25 -td 0.25. The hg38.UCSC.add_miR.160920.refgene.mat file is included with GISTIC2.\nThe threshold for the residual q-value for both broad and focal amplifications and deletions was set at 0.10. To identify driver events, we focused on genes within significant amplifications and deletions <= 10Mb in size. To remove artefactual recurrent SCNAs, we further removed GISTIC2 wide peaks with an overlap of more than 40% with regions known to be problematic for sequencing and/or read mapping, as defined by the Genome in a Bottle Consortium benchmark union (v3.3; https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/genome-stratifications/v3.3/GRCh38@all/Union/GRCh38_alldifficultregions.bed.gz). Finally, we also excluded a wide peak if the list of samples predicted to have a particular SCNA in GISTIC2 did not concur with the original ASCAT results.\ni.e. to be defined as a significant focal copy number event a region needs to have a GISTIC2 residual q-value < 0.1, a peak size > 100,000bp, show concordance with the ASCAT CN call of >= 0.75, and for the region to not overlap with a Genome In A Bottle \"difficult regions\" by >40%.\nThe manuscript was amended to clarify the above (please see revised Manuscript).\nCompelling examples of the raw copy number data for example for HRNR and KIT amplifications should be shown for this section to be more compelling, ideally with additional supporting information nominating these genes (FISH, IHC).\nLooking at absolute CN calls for Hornerin (HRNR) and KIT, we can clearly see these 2 genes are amplified as shown below:\nHRNR\nPD42480a\nPD42532a\nPD42542a\nKIT\nPD47219a\nFigure R4 - Absolute CN calls of KIT and HNRN in sebaceous tumours.\nWe do not have the material to do FISH and based on the RNA-Seq data we know that these genes are expressed so we don’t feel that IHC will be informative and at best such an analysis would be semi-quantative.\nA better integration with mutation data is desirable when nominating specific genes based on copy number change. E.g. how would an amplification of GNAS be plausible in this tumor type, in which no activating mutations of GNAS are seemingly present?\nWe looked at the GNAS amplification (as well as the other candidate oncogenes on the same cytoband). None of the amplified GNAS sample had a co-occurrent somatic mutation of GNAS (same for SALL4, NFATC2 and PTK6). Please see “NatureComs_GNAS” in “SupplementaryMaterial” folder.\nSeveral cancer-related genes are frequently amplified but not commonly mutated in tumours. Amplification generally refers to the increase in the number of copies of a particular gene, which can lead to overexpression/mis-expression of the encoded protein and in this way contribute to tumorigenesis. This does not have to be with concomitant mutation of the gene.\nHere are some key examples:\n1. ERBB2 (HER2/neu)\n- Function: Encodes a receptor tyrosine kinase involved in cell growth and differentiation.\n- Amplification: Commonly amplified in breast cancer, gastric cancer, and other solid tumours.\n- Clinical Significance: Overexpression of HER2 is associated with aggressive tumour behaviour and poor prognosis. HER2-targeted therapies (e.g., trastuzumab) are effective in treating cancers with HER2 amplification.\n2. MYC (c-MYC)\n- Function: Encodes a transcription factor that regulates genes involved in cell proliferation, growth, and metabolism.\n- Amplification: Often amplified in various cancers, including breast cancer, lung cancer, and neuroblastoma.\n- Clinical Significance: MYC amplification is associated with high tumour grade, aggressive behaviour, and poor prognosis.\n3. MDM2\n- Function: Encodes a protein that inhibits the tumour suppressor p53.\n- Amplification: Frequently amplified in sarcomas, glioblastomas, and certain carcinomas.\n- Clinical Significance: MDM2 amplification leads to p53 inactivation, contributing to tumour progression. It is a potential target for cancer therapies.\n4. CCND1 (Cyclin D1)\n- Function: Encodes a protein that regulates cell cycle progression.\n- Amplification: Commonly amplified in breast cancer, head and neck cancers, and lymphomas.\n- Clinical Significance: Cyclin D1 overexpression due to gene amplification drives cell proliferation, making it a target for therapeutic intervention.\n5. EGFR (Epidermal Growth Factor Receptor)\n- Function: Encodes a receptor tyrosine kinase involved in cell growth and survival.\n- Amplification: Frequently amplified in glioblastomas, non-small cell lung cancer, and other tumours.\n- Clinical Significance: EGFR amplification is associated with increased cell proliferation and survival, and EGFR-targeted therapies are used in treating certain cancers.\n6. MET\n- Function: Encodes a receptor tyrosine kinase involved in cellular growth, motility, and survival.\n- Amplification: Often amplified in gastric cancer, lung cancer, and papillary renal cell carcinoma.\n- Clinical Significance: MET amplification can lead to constitutive activation of the receptor, driving oncogenesis. Targeted therapies against MET are being developed.\n7. CDK4\n- Function: Encodes a cyclin-dependent kinase that plays a key role in cell cycle progression.\n- Amplification: Commonly amplified in melanoma, sarcoma, and glioblastoma.\n- Clinical Significance: CDK4 amplification promotes cell cycle progression and is a target for CDK4/6 inhibitors in cancer therapy.\n8. FGFR1 (Fibroblast Growth Factor Receptor 1)\n- Function: Encodes a receptor tyrosine kinase involved in cell proliferation, differentiation, and angiogenesis.\n- Amplification: Amplified in lung squamous cell carcinoma, breast cancer, and other tumours.\n- Clinical Significance: FGFR1 amplification drives tumorigenesis, and FGFR inhibitors are being explored in clinical trials.\nThese genes are frequently amplified in various cancers, only occasionally with mutation, and their amplification often plays a critical role in driving tumour progression, making them key targets for cancer therapies.\nThe wording and grammar in these paragraphs is awkward e.g. ‘In contrast, SC-Es showed significantly chromosomes 1 and 20 amplifications, and significantly 17p deletions, which includes TP53 (7/36 tumours, of which 6 also had a somatic TP53 mutations’\nThe entire edited paragraph is below: “At a chromosome arm level, SC-Os showed the greatest number of significant recurrent SCNAs, including chromosomes 1q, 8q, 19p and 20 amplifications, as well as chromosome 4 and 13q deletions (Supplementary Figure 10A-D), with the 13q deletion encompassing RB1 and BRCA2 (8/15; one of them being somatically mutated for RB1) (Supplementary Figure 10E). In contrast, SC-Es showed significant chromosomes 1 and 20 amplifications, and significant 17p deletions, which encompass TP53 (7/36 tumours, of which 6 also had a somatic TP53 mutations; Supplementary Table 11).”\n(please see revised Manuscript).\n6) Additional consideration and discussion should be given to the PAK2 fusions. The fusions with PAK2 at the 5’ end appear to be novel/unprecedented. Is there a reciprocal fusion product that was perhaps missed?\nNotably, PAK2 fusions (PAK2::DLG1, PAK2::LRIG1 and DLG1::PAK2) were identified Supplementary Figure 11A and Supplementary Table 14) with the reciprocal event DLG1::PAK2 observed in the case with a PAK2::DLG1 fusion. We did not identify a reciprocal event in the tumour with the PAK2::LRIG1 fusion which might suggest it is expressed below the level of detection or it is not generated. Importantly, our analysis suggests that PAK2 fusions are associated with benign tumours and with sebaceous differentiation as was recently reported in apocrine poroma as well as in two SMs and one SA, without specific features, in our cohort (Supplementary Table 14).\n7) The finding of an identical fusion of RCOR1::GRHL2 needs additional discussion as it implies a clonal relationship between the two SAs of this patient.\nOut of the two sebaceomas samples from this patient sample PD43859a did not passed DNA QC. Thus, to explore the question of a shared clonal origin further we compared the mutations found in the DNA/exome sequence of PD43859c to variants called from the transcriptome sequence of PD43859a.\nThis analysis revealed:\n Of the 7 sites called as somatic variants from PD43859c with sufficient coverage in PR43859a 5 / 7 variant calls were concordant. i.e. they appear as somatic variants in both PR43859a and PD43859c.\n If we consider all sites, at least 8 / 15 ((5+3) / 15) variants in ‘c’ are found in ‘a’; however, 5 sites had coverage < 10X and therefore may be false negatives.\nCollectively, these observations do suggest a shared origin for a and c. However, these two samples were both benign and collected with a time gap of 3 years between each excision. This time gap means that it is unlikely to be a recurrence. We therefore wonder if these tumours have a shared, possibly post-zygotic, origin with this specific patient noted as having multiple (many) sebaceous cysts in their clinical notes. This is an intriguing observation indeed.\n8) The finding on CIN primarily in pMMR tumors is interesting and mirrors the findings in CRC. This indicates, as you state, that genomic instability is necessary for the malignant transformation. Any insights on why and how?\nComparisons of CIN (measured in terms of numbers of gains and losses of chromosome arms or long DNA segments therein) in pMMR and dMMR colorectal cancers (CRC) has shown a clear difference with pMMR CRCs showing significantly more CIN than dMMR CRCs. These data are consistent with the hypothesis that at least one form of genomic instability, either CIN as seen in the majority of CRC that are pMMR, or dMMR is required to allow transformation to malignancy. Here, we hypothesise similar mechanisms occurring in sebaceous tumours with CIN in pMMR or MSI in dMMR cases. Underlying CIN there can be increased replication stress and defects in mitotic spindle assembly, errors in chromosome cohesion and segregation, centrosome amplification, and/or dysfunctional checkpoints. Telomere dysfunction and epigenetic mechanisms may also contribute.\nReferences:\nDNA Mismatch Repair Deficiency in Sporadic Colorectal Cancer and Lynch Syndrome.\nPoulogiannis G, Frayling IM, Arends MJ. Histopathology 2010; 56: 167 - 179.\nPrognostic relevance of DNA copy number changes in colorectal cancer.\nPoulogiannis G, Ichimura K, Hamoudi R, Luo F, Leung SY, Yuen ST, Harrison DJ, Wyllie AH,Arends MJ. Journal of Pathology 2010; 220: 338 - 347.\nLynch Syndrome – Cancer Pathways, Heterogeneity and Immune Escape\nSeth S, Ager A, Arends MJ, Frayling IM.\nJ Pathology 2018 246(2):129-133. doi: 10.1002/path.5139. PMID: 30027543\nMolecular Pathology of Lynch Syndrome.\nCerretelli G, Ager A, Arends MJ, Frayling IM.\nJ Pathol. 2020 Apr;250(5):518-531. doi: 10.1002/path.5422. PMID: 32141610\nLynch Syndrome. Frankel W, Arends MJ, Frayling IM, Nagtegaal ID. In: “World Health Organization Classification of Tumours of the Digestive System”. 5th Edition. IARC Press: Lyon, France, 2019. Chapter 14, pp515-521. Genetic Tumour Syndromes of the Digestive System. ISBN 978-92-832-4499-8.\nGenetic tumour syndromes associated with skin malignancies: Muir-Torre syndrome.\nTetzlaff M, Arends MJ, Wood BA. 2024, World Health Organization Classification of Skin Tumours. 5th ed. IARC Press, Lyon, France. Chapter 9, Section 9.7 (in press).”\nWe thank Reviewer 1 for their comments.\nReviewer #2 (Remarks to the Author):\nGeneral comments: expertise in multi-omics integration\nIn this multi-center international study, Ferreira et al. present the largest exome and transcriptome analysis of sebaceous tumors (ST) to date. This allowed investigators to identify mutated genes under positive selection (drivers) for the first time in this tumor type, which provides an incredibly valuable resource to the cancer research community. In my opinion, this is an extremely important study for the ST and cancer genomics fields.\nWe thank Reviewer 2 for supporting our project and highlighting the fact this is an incredibly valuable resource and extremely important study.\nThe initial data processing, including variant calling, mutation signature extraction, and copy number profiling, was conducted using well-established software and pipelines. However, several downstream analyses are problematic, and the presentation of the results can be significantly improved to enhance clarity and better guide the reader. Furthermore, in my opinion, parts of the neoantigen and ligand-receptor analyses are too speculative, and they lack additional data or experimentation for validation. These sections should be modified for clarity and rigor, and their results presented more cautiously or the authors may want to consider omitting these analyses from the study.\nWe thank Reviewer2 for sharing their concern about the neoantigen and ligand-receptor analyses. We agree these analyses are purely speculative and that follow-up functional analysis are required to confirm them. We have significantly reduced or removed these analyses or moved them to the supplementary information (Figure 7 has become Supplementary Figure 15, for example) and we have rephrased the text to make clear follow-up studies will be required for validation.\nMy recommendation is to improve the whole exome and transcriptome analysis as outlined in detailed below and re-consider inclusion of the ligand-receptor analysis and parts of the neoantigen analysis.\nWe thank the reviewer for their comments which we have addressed in full below.\nConcerns/suggested revisions:\nSomatic landscape of sebaceous tumours\nDriver Identification Analysis:\n[1] The paper reports significantly mutated genes based on a single software tool (dNdScv), without the necessary validation steps (e.g. filtering for genes that are not expressed or have low VAF), raising concerns about false positives that are ubiquitous in the context of high TMB tumors. (For example, SPANXB1 is reported as a driver in Figure 2, but most mutations in this gene appear to have a suspiciously low VAF). To ensure robustness, the authors should filter their results to exclude likely false positives and, importantly, provide additional supporting evidence for the novel candidate drivers:\n● What is the variant allele frequency (VAF) distribution for each gene?\nAs requested, we have made a plot showing the VAF distribution for each driver gene (please see Figure R5). We can observe that SPANXB1 is not outside the distribution, while ZNF750 (VAF (mean): 0.479661) is the only statistically significant outlier based on the interquartile range (IQR) method (please see Figure_R5 in “SupplementaryMaterial” folder).\ni.e. only ZNF750 is outside the distribution with the VAFs for this gene being higher than for other drivers.\nFigure R5 - Driver genes and variant allele frequencies (VAF) in sebaceous tumours. (A) Barplot showing the distribution of the VAF across the driver genes. Each dot is a tumour.\n● Do these variants mostly originate from a specific signature or positional site?\nAs requested, we have made a plot showing the VEP variant classifications across the 28 driver genes identified in this study (Figure R6/Supplementary Figure 4A). Further, we also provide a figure illustrating the profile of nucleotide changes found in these genes (Figure R7 see below). This analysis revealed that most driver genes are characterised by a range of variant types, except HRAS, TMEM86A and SGCB, which contain only missense mutations suggesting they are oncogenes.\nFigure R6/Supplementary Figure 4A - Driver genes and associated variant types in sebaceous tumours. Barplot showing the proportion of the different variant types across the driver genes.\nFigure R7/Driver genes and associated nucleotide changes in sebaceous tumours. Barplot showing the proportion of the different nucleotide changes across the driver genes.\n● What proportion of variants are seen in the RNA-Seq data for each gene?\nAs suggested, we checked for the presence of the somatic variants call from our exome data in the RNA-seq data generated on each case. We can observe that most of the somatic mutation found in the DNA are present/expressed in the RNA-seq.\nVariant Mutations with Total mutations Mutations Variants with <=\ntype coverage > 10x in observed found 10x coverage\nRNAseq in RNAseq in RNAseq in RNAseq\nSNV 676 667 98.7% 202\nDNV 15 15 100% 5\nDEL 414 407 98.3% 104\nINS 128 133 96.2% 55\n● Are there significant co-occurrences of loss-of-function mutations with loss of heterozygosity (LOH) for the putative tumor suppressor genes (TSGs)?\nFor each driver genes, we have calculated the percentage of LOH and made a new plot for(below):\nFigure R8 - Supplementary Figure 4C - Driver genes and associated copy number alterations in sebaceous tumours. Barplot showing the proportion of the different copy number alterations across the driver genes.\nWe can see that some of the established oncogenes we detected were amplified (such as KRAS), while loss and cnLOH of TSGs (such as MSH2 and RB1) was identified. Surprisingly, SPANXB1 is frequently affected by cnLOH.\n[1.1] The authors of the dNdScv software have excluded hypermutator samples (>500 mutations) from their main analysis (iMartincorena et al, 2017, PMID:29056346). Such precautions are not mentioned in the method section of the current manuscript. Did they authors try running dNdScv exclusively on tumors with fewer than, say, 1000 or 2000 mutations to see how it affects the results?\nWe thank Reviewer2 for this comment. We consulted the lead-author of the dNdScv software, Inigo Martincorena, who is a member of our research division. He explained to us that “500 mutations” was an arbitrary number randomly chosen without any statistical testing/justification. As noted, we performed the dNdScv analysis on all mutations but in response to the reviewers comment we also provide the analysis performed with other cut-offs. More specifically, we ran dNdScv removing samples with >2000, >1000 and >500 mutations. Critically, regardless of the cut-off applied, 88%, 60% and 40%, of driver genes were retained when running the analysis across the entire cohort with driver gene detection appearing to mainly be affected by the number of samples used for the analysis at each TMB level. By providing this analysis in the revised manuscript readers can explore the data across these threshold levels. See “Supplementary_Table_5”).\n[1.2] There are 24 samples with no matched normal. It is important to exclude these samples from the dNdScv analysis as they could contain germline variants misclassified as somatic (I’ve noticed VAF modes near 0.5 for many of these samples). Could the authors specify which samples have been included in their analyses, along with any exclusion criteria?\nWe thank Reviewer2 for this comment. The additional 24 samples were only included in the plots in Supplementary Figure 3. No other analyses included them.\n[1.3] Possible additional figures for the driver analysis (in main or supplementary): A barplot showing the relative proportion of each variant class (missense, frameshift stop gain etc) in each gene or selected genes. A scatter plot showing dN/dS:missense vs dN/dS:Lof similar to Figure 2B in Martincorena et al.(2017, PMID: 29056346). These figures would help support the analysis results and importantly inform about the action mode of the candidate drivers (i.e. oncogene vs.TSG).\nAs requested, the two new plots were generated and included in Supplementary Figures 4A-B in the revised manuscript.\nFigure R9/Supplementary Figure 4A - Driver genes and associated variant types in sebaceous tumours. (A) Barplot showing the proportion of the different variant types across the driver genes.\nFigure R10/Supplementary Figure 4B - Driver genes and associated variant types in sebaceous tumours. (B) Scatterplot showing missense versus truncating mutation profile across the driver genes.\n[1.4] It would be helpful if the authors cross-reference their gene list with the results of some pan-cancer analyses, such as Bailey et al. (2018, PMID: 29625053) or Martincorena et al. (2017, PMID: 29056346). This information could be easily summarized in a supplementary figure having genes as rows and cancer types as columns, with each tile colored according to their report status (e.g. oncogene, potential oncogene, tsg and potential tsg, which is available in Supplementary Table 1 of Bailey et al.). It would be useful for the readers for genes that haven’t been reported as significantly mutated, what does the literature say about them that could support their role as an oncogene or tumor suppressor?\nWe made the plot requested cross-referencing our results with those of Bailey et al. (2018, PMID: 29625053) and we have added it to Figure R11/Supplementary Figure 4D.\nFigure R11/Supplementary Figure 4D – Comparison of mutation consequences in sebaceous tumours driver genes and across TCGA cancers. Tileplot showing the alteration of the driver genes across sebaceous tumours and pan-cancer161. The function of each gene has been defined previously161. Abbreviations: NA, not defined; onco, oncogene; onco.tsg, oncogene and tumour suppressor genes; tsg, tumour suppressor gene.\n[1.5] I speculate the authors may have missed the identification of several drivers by testing all genes as they did impose a burden on significance levels after multiple testing corrections. In addition to their global analysis, I think it may be worth performing a restricted hypothesis testing (RHT), such as in Lawrence et al. (2014, PMID: 24390350) by considering only previously reported genes, for example those from Bailey et al. (2018), and those reported as drivers in sebaceous tumors. This additional analysis would be straightforward, requiring only the subsetting of the global results to these genes before applying multiple testing corrections. This may be valuable to the readers.\nAs requested, we have performed restricted hypothesis testing (RHT) with the dNdScv algorithm using the list of established driver genes reported by Bailey et al., (please see “Supplementary Table 5”). This analysis increased the detections of candidate driver genes as predicted (19 new genes were identified as driver while the previous reported driver genes were kept). We have drawn another oncoplot with this new list of driver genes and this can be found in Figure R12/Supplementary Figure 5.\nFigure R12/Supplementary Figure 5 - Somatic mutational landscape (driver genes post restricted hypothesis testing (RHT)) of sebaceous tumours across tumour subtypes. Oncoplot showing somatic mutations in the driver genes after RHT (ranking based on the percentage of samples altered (decreasing order)) of sebaceous tumours (STs). This includes the 44 driver genes identified by dNdScv when taking all four ST subtypes together, and performing RHT. Patterns of mutation of the top 4 driver genes (NOTCH1, RREB1, TP53 and HRAS) allows division of the STs into five molecular subtypes (A-E) across the four ST subtypes. ‘Site’ is the body location of the tumour. The ‘PD’ number is the sample number, with those in red font being from patients that developed metastases. TMB is in mutations/MB, including SNVs, MNVs and indels). Abbreviations: NV, not available.\n[1.6] The main oncoplot could be enhanced in several ways.\n• Adding copy number track for the most frequently affected genes (deleted TSG or amplified oncogenes).\n• Some drivers with less supporting evidence or of lower interest could be removed from the main figure for simplicity and clarity (but shown in a full supplementary oncoplot).\n• Annotation tracks with tumor purity, RNA subtypes and mutational signature groups could be added.\n• Mode of action of each driver (TSG vs oncogene) could be highlighted.\nAs requested, we have added a track highlighting the proposed mode of action of each driver in Figure R13/Figure 2 (red for oncogene and blue for TSG). However, after trying to add all the requested information, the figure became unreadable, so we elected to retain this information in other figures and in the supplementary Tables. Regarding the SCNAs that overlap with the driver genes, the only gene significantly recurrently deleted is RB1, which can be seen in Figure 5 as a significant peak for SC-E (cytoband 13q14.3) and as a broad peak/arm in SC-O in Supplementary Figure 10 (Arm 13q). Therefore, we did not include this track in the figure.\nFigure R13/Figure 2 – The somatic mutational landscape of sebaceous tumours across subtypes.\n(A) Oncoplot showing somatic mutations in the driver genes of sebaceous tumours (STs). This includes the 25 driver genes identified when taking all four ST subtypes together, and three additional driver genes when considering only specific ST subtypes. Gene names in red denote oncogenes, genes in blue denote tumour suppressors, while genes names in both colours are those reported to have both functions depending on the text. Patterns of mutations in the top 4 driver genes (NOTCH1, RREB1, TP53 and HRAS) allows division of the STs into five molecular subtypes (A-E) across the four ST subtypes. ‘Site’ is the body location of the tumour. The ‘PD’ number is the sample number, with those in red font being from patients that developed metastases. (B) Mutual exclusive and co-occurring interactions based on DISCOVER. (C) Comparison of the tumour mutation burden (TMB) of the sebaceous tumours (STs) against 33 TCGA cohorts (TCGA tumour types are shown using standard abbreviation). The ST cohort are shown all together (blue font), or separated by subtype and mismatch repair (MMR) status (green and red font). Abbreviations: MMR, mismatch repair; dMMR, deficient in MMR; pMMR, proficient in MMR; NV, not available; TMB, tumour mutational burden (in mutations/Mb, including SNVs, MNVs and indels).\n[1.7] The authors only report the gene names with their corrected p-values in the supplementary material. It is important to include the complete dNdScv output in the supplementary material to provide a comprehensive view of the analysis and facilitate further validation and reproducibility of results.\nAs requested, we provided the data into “NatureComs_dNdScv_Output” in “SupplementaryMaterial” folder) and we will also make this available via the GitHub repository.\n[1.8] On lines 199-200, the authors write, \"Interestingly, when considering the ST subtypes individually, HRAS was only identified as a driver gene in the benign tumors.\" What is the uncorrected p-value for HRAS in the SC-E and SC-O subgroup analyses? Would it pass RHT? Even if not significant, the oncoplot shows several HRAS mutations in SC-E and SC-O, most of which are known activating mutations at G12/13 and Q61. Moreover, HRAS harbors no silent mutations in these tumors, strongly suggesting it is an actual driver, albeit less frequently mutated. Accordingly, it would be more appropriate to drop this statement or at least revise it to \"HRAS was not reported as significant in (...) probably due to its lower mutation frequency and the smaller size of these groups (...),\" and mention the presence of these hotspot mutations. In general, the absence of a significant result for a particular gene in a subgroup analysis cannot be used on its own to suggest an absence of “driver” activity.\nWe thank Reviewer2 for this comment and have amended the text as follow:\n“Interestingly, when considering the ST subtypes individually, HRAS was only identified as a statistically significant driver gene in the benign tumours, although hotspot HRAS mutations were also observed in malignant subtypes, but these mutations occurred at a lower mutation frequency in these smaller cohorts (Supplementary Table 5).”\n[1.9] In the discussion, lines 491-492, the authors state: “Use of somatic mutational profiling may help in the differential diagnosis of benign versus malignant STs, which can be challenging.” It would be helpful if more clinical context could be provided here, particularly around the standard histopathological diagnostic workflow. What are the gold-standard IHC stains, histopathological correlates, or clinical decision algorithms currently used to differentiate between the four types of sebaceous tumors? A deeper discussion on how mutational testing could provide clinical input would be ideal, given the clinical correlation angle in the discussion. Figure 2 does not clearly identify a mutational profile with high sensitivity or specificity for diagnosis. Could the authors consider using classifiers (e.g., Random Forests) and measure their performance with cross-validation?\nWe thank Reviewer2 for highlighting the importance of correlating clinicopathological and genetic aspect. Our team includes Professor Thomas Brenn who is a world-leading dermatopathologist and author of several textbook such the WHO blue book for skin and the McKee – Pathology of the skin, and a noted expert in adnexal tumours. At present, the gold standard to make the diagnosis of these sebaceous tumours is based on histology only, there is no core set of tests/IHC stains. Through our paper, we suggested markers/tools to better detect Lynch syndrome thanks to clinical and histopathological variables (please see “Figure 1”). The use of somatic mutation profiling may help in the diagnosis of benign and malignant STs. In the first paragraph of the discussion we say the following:\n“Several findings in this study may prove clinically relevant in the diagnosis of ST and LS. Firstly, somatic mutational profiling may help in the differential diagnosis of benign versus malignant STs, which can be challenging. For example, the malignant STs were driven by chromosomal instability and typically showed biallelic TP53 inactivation, as well as TP53 mutations co-occurring with ZNF750 and/or RB1 mutations. In addition, dMMR SC-E was characterised by mutations in NOTCH1, RREB1, and DRD5, whereas SC-O was characterised by significant focal amplifications encompassing GNAS/SALL4 and LYN/PLAG1, significant broad 13q deletions (encompassing BRCA2 and RB1), and 8q amplifications (encompassing MYC), with amplification of the MYC proto-oncogene previously reported by targeted sequencing and FISH 11. Lastly, both SC-E and SC-O showed significant amplification of chromosomes 1q21.3 and 20. The former amplification encompasses HNRN (Hornerin), which is part of the epidermal differentiation complex on 1q21. Hornerin is expressed by the ocular surface (corneal and conjunctival aspects, including meibomian and Zeis glands), the lacrimal apparatus and the skin, and has antimicrobial and protective (keratinisation) functions 47,48, reinforcing its potential role in sebaceous carcinoma. STs, with the exception of SC-O, can be the first sentinel sign of LS and 3 clinical criteria and 4 histopathological criteria – including (i) the presence of multiple STs, (ii) STs being located outside the head and neck area, (iii) the presence of a LS-associated cancer, in the clinical model, and (i) the absence of epidermal connection, (ii) the delimitation of the ST by a collarette, (iii) a KA-like architecture, and (iv) the diagnosis of SM in the histopathological model, together with loss of MMR expression by IHC – were found to significantly associate with a LS diagnosis (the clinical model being a better predictor) and thus should raise the index of suspicion that could provoke a referral for germline genetic testing. Of note, while the idea of generating a classifier is one we feel is very validate, a comparable dataset of STs to the one we describe in our paper does not exist which makes validating such a tool problematic and thus we feel this analysis is premature. Further, this study highlights the importance of recording the precise location of SC-O tumours, and specifically restricting use of the term SC-O to those lesions arising on the eyelid. Indeed, the sample coming from the canthus showed a different phenotype from the other SC-O and was characterised by a UV mutational signature, rather than being pauci-mutated, the pattern observed in SC-Os.”\nMutually exclusive and co-occurring mutations:\n[2] The “mutual exclusions and co-occurring interactions” analysis is problematic as it combines the results of a Fisher’s test with the output of DISCOVER, a tool specifically developed to replace the inappropriate usage of the Fisher’s test in such analyses. The Fisher’s test does not control for TMB, leading to a high number of false positive co-occurring interactions (as shown in Figure 2B). In contrast, DISCOVER accounts for TMB in each sample and gene, thus limiting these spurious results. Combining the results from both approaches does not make sense, as exemplified by the contradictions between the interactions reported in Figure 2B and the DISCOVER analysis output provided in Supplementary Table 6 (e.g., significant co-occurrence between ZNF750 and BICRA while Supplementary Table 6 reports significant mutual exclusion). I would suggest that the authors should rely solely on DISCOVER (or similar approaches controlling for TMB), and completely exclude the Fisher’s test from this analysis.\nWe thank Reviewer2 for sharing their concerns about the use of Fisher test. As suggested, we have removed this analysis (please see Figure R13/Figure 2). We also ran DIALECT (Shuaibi A et al. 2024; https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11071465/) on the cohort. The DIALECT software distinguishes between missense and nonsense mutation by design and claims to be superior to DISCOVER. However, even after talking to the authors and re-writing some of their code we were unable to deploy this algorithm.\nFigure R13/Figure 2 – The somatic mutational landscape of sebaceous tumours across subtypes.\n(A) Oncoplot showing somatic mutations in the driver genes of sebaceous tumours (STs). This includes the 25 driver genes identified when taking all four ST subtypes together, and three additional driver genes when considering only specific ST subtypes. Gene names in red denote oncogenes, genes in blue denote tumour suppressors, while genes names in both colours are those reported to have both functions depending on the text.\nPatterns of mutations in the top 4 driver genes (NOTCH1, RREB1, TP53 and HRAS) allows division of the STs into five molecular subtypes (A-E) across the four ST subtypes. ‘Site’ is the body location of the tumour. The ‘PD’ number is the sample number, with those in red font being from patients that developed metastases. (B) Mutual exclusive and co-occurring interactions based on DISCOVER analysis. (C) Comparison of the tumour mutation burden (TMB) of the sebaceous tumours (STs) against 33 TCGA cohorts (TCGA tumour types are shown using standard abbreviation). The ST cohort are shown all together (blue font), or separated by subtype and mismatch repair (MMR) status (green and red font). Abbreviations: MMR, mismatch repair; dMMR, deficient in MMR; pMMR, proficient in MMR; NV, not available; TMB, tumour mutational burden (in mutations/Mb, including SNVs, MNVs and indels).\nMutational signatures of sebaceous tumours:\n[3] The mutational signature analysis was not comprehensively explored and presented, and a global view of the cohort is lacking, specifically related to the relative contribution of the main signatures in each sample. The presentation of the results in Figure S5 show too many signatures simultaneously, and additional figures and analyses would be necessary to support the claims and observations made in the main text and discussion. In my opinion, an improved global view of the mutational signatures would warrant inclusion as a main figure.\nWe thank Reviewer 2 for their comment. We have log-transformed the data which makes the SBS signatures easier for the reader to see. This figure is available in Figure R14/Figure 4A of the revised manuscript as suggested by the reviewer.\nFigure R14/Figure 4A – Mutational signatures of sebaceous tumours. Bar plot showing the mutational signature profiles, including SBS (single-base substitutions), ID (insertions and deletions) and DBS (doublet-base substitutions) signatures, across each of the four sebaceous tumour subtypes. ‘PD’ numbers refer to the DNA sample identifiers while the asterisk denote samples with/without SBS or ID signatures. A black asterisk or a PD number in black (*/PD A) indicates that a signature was detected in that sample, while blue (*/PD NV) a or red asterisk or PD number (*/PD NA) indicates samples where no mutational signature was detected or samples were not analysed, respectively. COSMIC signature groups: SBS: clock-like signatures (SBS1 and SBS5), APOBEC signatures (SBS2 and SBS13), tobacco signature (SBS4). UV signatures (SBS7a, SBS7b, SBS7d), POLE signatures (SBS10a and SBS10b), dMMR signatures (SBS15, SBS20, SBS21 and SBS26; SBS20 is linked to dMMR and POLD1 mutation), azathioprine signature (SBS32), MUTYH signature (SBS36), unknown signature (SBS28). ID: slippage signatures (ID1 and ID2), dMMR signature (ID7). DBS: UV signature (DBS1).\n[3.1] Figure S5 needs clarification. The graphical legend mentions */PD A, */PD NA, and */PD NV in three different colors, but none of these are explained in the textual legend. There is a reference to the blue and red labeling of the sample names, but this is confusing since \"blue\" samples (no significant signature) actually show some signatures. Additionally, the legend states that \"one in red text was not included,\" but half of the samples are in red.\nAs suggested by Reviewer 2 we have clarified these symbols (*/PD A, */PD NA, and */PD NV) in the revised manuscript.\nFigure 4A – Mutational signatures of sebaceous tumours. Bar plot showing the mutational signature profiles, including SBS (single-base substitutions), ID (insertions and deletions) and DBS (doublet-base substitutions) signatures, across each of the four sebaceous tumour subtypes. ‘PD’ numbers refer to the DNA sample identifiers while the asterisk denote samples with/without SBS or ID signatures. A black asterisk or a PD number in black (*/PD A) indicates that a signature was detected in that sample, while blue (*/PD NV) or red asterisk or PD number (*/PD NA) indicates samples where no mutational signature was detected or samples that were not analysed, respectively. COSMIC signature groups: SBS: clock-like signatures (SBS1 and SBS5), APOBEC signatures (SBS2 and SBS13), tobacco signature (SBS4). UV signatures (SBS7a, SBS7b, SBS7d), POLE signatures (SBS10a and SBS10b), dMMR signatures (SBS15, SBS20, SBS21 and SBS26; SBS20 is linked to dMMR and POLD1 mutation), azathioprine signature (SBS32), MUTYH signature (SBS36), unknown signature (SBS28). ID: slippage signatures (ID1 and ID2), dMMR signature (ID7). DBS: UV signature (DBS1).\n[3.2] The method section states that SigProfilerExtractor was used for de novo extraction, considering solutions with 1 to 10 signatures by NMF. However, Supplementary Figure 5 shows 21 signatures (17 SBS, 3 ID, and 1 DBS), and Supplementary Table 8 reports mutation counts for 27 signatures. This discrepancy needs clarification. It is essential to explain the criteria for determining the final set of signatures and how they relate to the initial NMF analysis to ensure consistency and transparency.\nWe thank Reviwer2 for their comment. For de novo extraction, we ran the analysis using 1 to 10 signatures. However, the next step is signature decomposition and signature assignment. At this step, mutations were assigned to 27 of the 257 signatures in the COSMIC catalogue. SigProfilerExtracter was run in exome mode using GRCh38 as the reference and opportunity genome and 500 replicates for non-negative matrix factorisation (NMF) with 1 to 10 signatures. The solution implemented for downstream analysis was the optimal solution provided by SigProfillerExtractor (referred to as the ‘suggested solution’). Results from signature decomposition and assignment were taken forward if the cosine similarity between the de novo extracted signatures and signatures reconstructed from assigned COSMIC signatures was >= 0.9.\nThe de novo signatures are decomposed and assigned to the known COSMIC signatures, which is why there are more than 10 COSMIC signatures.\n[3.3] I strongly suggest the authors include a main figure summarizing the mutational signature analysis. For example, a clustered heatmap of each signature log10 counts (as rows) in each sample (columns) clearly shows the 3 main mutational profiles within this cohort: MMR (SBS15 + ID2 + ID7 + high SBS1 + high SBS5), UV (DBS1 + SBS7a + SBS7b + low SBS1 + low SBS5) and others (low SBS1 + low SBS5). Annotation tracks for samples could include the MMR status, log10TMB, diagnosis and presence/absence of LoF mutations in MSH2 & MSH6. An annotation track for the signatures could inform the reader about their aetiology (UV, MMR etc). Such heatmap can be easily produced with R’s pheatmap package and the ward.D2 clustering method. To make it more concise and suitable for a main figure, some highly correlated or same aetiology signatures could be merged (e.g. SBS7a & SBS7b) and low contributing signatures could be dropped and deferred to supplementary figures (if necessary). Alternatively, the authors could use the W matrix from NMF or apply PCA on the 96-SBS mutation proportion matrix to create a scatter plot of all samples colored by cancer types, MMR status and other variables of interest.\nWe thank Reviewer2 for their suggestion and have drawn the suggested heatmap and added it to Figure R15/Figure 4B.\nFigure R15/Figure 4B - Clustered heatmap of mutational signatures (log10 counts) in sebaceous tumours. Heatmap showing the mutational signature which tend to cluster by aetiology (left border of the plot). Colour gradient going from high mutations count (in red) to low mutational count (in yellow), to no mutational count (in light blue), to sample not included into the analysis (dark blue). COSMIC signature groups: SBS: clock-like signatures (SBS1 and SBS5), APOBEC signatures (SBS2 and SBS13), tobacco signature (SBS4). UV signatures (SBS7a, SBS7b, SBS7d), POLE signatures (SBS10a and SBS10b), dMMR signatures (SBS15, SBS20, SBS21 and SBS26; SBS20 is linked to dMMR and POLD1 mutation), azathioprine signature (SBS32), MUTYH signature (SBS36), unknown signature (SBS28). ID: slippage signatures (ID1 and ID2), dMMR signature (ID7). DBS: UV signature (DBS1).\n[3.4] Lines 268-269, the authors write, “As expected, in the dMMR cases, the high TMB was attributed to dMMR-associated mutational signatures.” What proportion of TMB is attributed to dMMR signatures? Could the authors please include a figure supporting this? It seems that other signatures (SBS1 & SBS5) contribute more to TMB than dMMR-associated signatures, even in dMMR samples.\nWe thank Reviewer2 for these suggestions and we made a figure supporting what we said (please see Figure R16/Supplementary Figure 8).\nFigure R16/Supplementary Figure 8 – Comparison between TMB, SBS signatures and ID signatures load in sebaceous tumours. Barplot showing the proportion of SBS and ID signatures related to TMB across the cohort.\n[3.5] Lines 163-166, the authors write, “Interestingly, indel rate rather than TMB, more accurately predicted MMR status (indel rate p-value < 2.2e-16 versus TMB rate p-value = 2.203e-08, Welch two samples t-test)”. I have a few concerns regarding this analysis. First, the considered data doesn’t meet this test's assumptions about normality. Second, p-values are not an appropriate measure of prediction accuracy. Third, according to the methods section, MMR status was in part determined by the indel rate (lines 772-774), which means using the indel rate to predict MMR status is circular reasoning. These issues should be addressed, for example by comparing ROC-AUCs and and an MMR classification that doesn’t rely on indel rate or TMB.\nWe thank Reviewer2 for their question. To test whether Indel rate or TMB are a better predictor of dMMR, we have used only samples that were defined as dMMR by having an IHC results suggesting loss of protein expression and also a biallelic alteration of a MMR gene (looking at somatic/germline mutations).\nWe then used a logistic regression for both Indel Rate and TMB and computed the ROC/AUC using leave-one-out cross validation (please see attached Figure R17/Supplementary Figure 1A).\nAs can be seen, IndelRate is a superior predictor over TMB:\nIndel Rate TMB\nMMR-p vs. MMR-d 0.9989177 0.8869048\nFigure R17/Supplementary Figure 2A - Receiver Operating Characteristic (ROC) curve for dMMR detection. Using logistic regression for both indel rate and TMB, and computing the ROC/area under the curve (AUC) using leave-one-out cross validation to define the best predictor for MMR status (pMMR versus dMMR), the AUC shows that indel rate (AUC = 0.999) is a better predictor than TMB (AUC = 0.887).\nThus, this figure shows indel are a better predictor than the TMB.\n[3.6] On lines 501-513, the authors write “this study highlights the importance of recording the precise location of SC-O tumours, and specifically restricting use of the term SC-O for only those arising on the eyelid. Indeed, the sample coming from the canthus showed a different phenotype from the other SC-O, being characterised by UV mutational signatures rather than being pauci-mutated”. Could the authors please include a figure and additional analyses or references supporting this claim, for example by showing differences in mutated drivers or gene expression that suggest it is more similar to SC-E than SC-O?\nWe thank Reviewer2 for this suggestion. We have made a new figure (please see Figure R18/Supplementary Figure 9) including SC-Es from the head and neck only, and the SC-Os, and the case from the canthus. We can easily see that the sample from the canthus shares mutations only found in SC-E such as CHD8, LATS1 and TMEM86A, as well as having a UV SBS mutational signatures.\nFigure R18/Supplementary Figure 9 - Mutational profile of SC-Es located on head and neck, and SC-Os. Comparison of SC-Es versus SC-Os across mutational signatures (A) and somatic mutations of drivers genes (B).\n[3.7] Lines 1377-1379: “SBS” should be “DBS”. Also, a typo in “mutationa”.\nWe thank Reviewer2 for pointing out these typos, which we have fixed (please see revised Manuscript).\n[3.8] In some other cancers, MMR status is linked to survival (e.g. Endometrial). Did the authors investigate whether in this cohort, MMR status is linked to survival in SC regardless of subtype? Authors note that sample size in SC was too small to measure but it may be helpful for readers to see this information anyways even if not statistically significant.\nWe thank Reviewer2 for this suggestion, however, only 4 patients in our cohort died from sebaceous tumours, of whom 2 metastatic SC-E, 1 metastatic SC-O and 1 SC-O with no reported metastasis where the patient died 4 months after diagnosis. Thus, we think a survival curve is therefore potentially misleading (please see Figure R19).\nFigure R19 – Survival curve in sebaceous tumours.\nIntegrative cluster analysis:\n[4] The results of this analysis are rather cryptic, and the methodological details raise concerns about their validity.\nWe agree that the integrative analysis was not properly described in the manuscript. In particular, we focused only on the clustering of tumours instead of the identification of molecular drivers, which is the main goal of this type of methods. We have now run two different methods: iCluster (https://doi.org/10.1093/bioinformatics/btp543) and MOFA (https://doi.org/10.15252/msb.20178124), both of them based on latent variable analysis and performed a comprehensive study on the robustness of the results. In the end, we decided to use MOFA, as the fitting process is more robust, and the results were easier to interpret. We must note that the goal of this method is not clustering, but identification of molecular drivers. We have included these data in the revised manuscript along with text and figures to explain the results.\nFirst, the expression data was combined with the mutation profiles of approximately 80 genes selected on the basis that they are frequently mutated in this cohort. Given that a typical tumor contains fewer than 10-15 driver mutations (referring here to Martincorena et al, 2017) and that these sebaceous tumors have huge TMBs (nearly half of them have a thousand mutations or more), we can be sure that most of the mutations in these 80 genes are passenger mutations. While a few relevant genes such as TP53 and MSH2 have been manually added and others like TTN have been manually removed, Supplementary Figure 4 suggests that the limited number of potentially meaningful mutations is still overwhelmed by the number of passengers.\nThe reviewer raises some interesting points. Indeed, the main goal of our integrative clustering is not so much to group tumours but to identify molecular driver events. To do that, both iCluster and MOFA methods identify sets of mutations and expression changes that co-occur in subsets of tumours. As both methods are computationally demanding, our first step in this analysis was feature selection. In the case of mutations, the reviewer is right that some of these tumours have a high mutational burden. We tried to identify a relatively small subset of genes that were potentially interesting, setting up filters based on the minimum number of samples that carry mutations. As the reviewer states, some of these could be passengers. We have added an additional filter testing for association of the mutation and the total mutational burden, removing those that were statistically significant. With this filter, we expected to remove those genes that were probably passengers. We also added manually some well-known driver genes. Of course, our final set of 84 genes might still contain some passenger genes, but both iCluster and MOFA will select only those that contribute to any driver event, as they contain regularisation to perform variable selection, removing uninteresting genes. In the case of expression, we used the 1,000 most expressed genes plus the genes present in the selected mutations.\nSecond, the method contains many steps with somehow arbitrary choices (e.g., a first selection step out of 16 fitted models followed by a second selection step for the number of k-means clusters), making it very difficult to assess the robustness of the findings. What happens if we use another fitted model (many had a similar balance of variance explained and a low number of clusters)? What happens if a different k is chosen or if we change the gene sets?\nThe optimal selection of the number of components is a well-known problem in Latent Variable Modelling and Cluster Analysis, and although several methods have been proposed, they often do not agree on the best solution. Starting with the selection of latent dimensions (that is, molecular drivers), MOFA is quite robust and in fact it recommends using the default number of 15. We checked the robustness of the models comparing the total variance explained and the minimum variance explained by any of the factors when we fitted a different number of latent factors in the model. As the plots below show, after 15 factors the increase in variance is minimal, therefore we used that value.\nWe also inspected the robustness in the model fitting using the recommended procedure in the vignette. We fitted our MOFA model using different random seeds and checked how similar the results were comparing the correlation of the factors obtained in each fit. We finally selected the best model.\nIt is also easy to quantify how much variance from each data type the latent factors explain, in order to focus on the most interesting drivers:\nIf we set up a threshold of minimum variance at 1%, we can retain for further study the first 10 latent factors. Indeed, they are the ones that contribute most to the tumours.\nThese drivers are characterised by both mutations:\nAnd gene expression:\nThese plots show the weights of each feature in each latent factor. The larger (or smaller) the x-value (weight), the more important the feature is in defining that latent factor. For example, having mutations in TP53 or ZNF750 will produce a high score in Factor 3. Looking at those weights we can also see that Factor 1 is explained mostly by mutations, and not that much by expression. In fact, that factor explains 23.8% of the mutation data and only 0.54% of the expression. It reflects high mutational burden, as all mutations have a similar positive score.\nThe following heatmap shows the weights of each molecular driver (latent factor, columns) in each tumour (rows):\nAs we can see each tumour presents a combination of these drivers, therefore it makes sense to try to cluster tumours based on their score. Using the top 10 factors we identified two potential solutions with 7 or 9 subgroups based on the silhouette index. These solutions were almost identical; the two new clusters were partition of the first and the third cluster. The following heatmap allows us to compare and choose the best partition. In fact, 9 clusters seem to identify better the subgroup of tumours that have high scores in the first factor (high mutational burden) and splits the third cluster into a group with high scores of factor number 4. We kept the solution with 9 integrative subgroups.\nWe describe now the 9 molecular subgroups:\nFor each integrative cluster (IntClust 1 to 9), the following data will be available:\n1) (n= ): the number of samples included into this IntClust;\n2) Subtype: which tumour subtype is mostly included into the IntClust;\n3) MMR: which MMR status represents the most the IntClust;\n4) Indel and SNV rate: which compare the mentioned IntClust against all the others;\n5) Proportion of cells: which compare the same cell type across the different Cluster and not the comparison between different cell type in a same Cluster;\n6) Proportion of mutated tumours: which proportion of samples in this IntClust are mutated for a specific gene;\n7) Pathway activation: activity of specific pathway in this IntClust;\n8) Copy number aberration: CNA profile of the IntClust.\n9) Mainly SC-O tumours.\nIntClust 1(n=13 tumours):\nA mutational-burden driven subgroup, mostly SM tumours with MMR deficiency and up-regulation of G2M checkpoint and E2F targets. High proportion of CD8+ T cells and fibroblasts.\nIntClust 2(n=15 tumours):\nIndel-driven cluster, but with lower indel levels. Mostly SM tumours. dMMR, MYC up-regulated and down-regulation in interferon response. High proportion of CD8+ T cells.\nIntClust 3(n=13 tumours):\nHigh indel burden, NOTCH1 and RREB1 mutations, dMMR but down-regulation of most hallmark pathways. High proportion of CD4+ T cells.\nIntClust 4(n=23 tumours):\nLower mutational burden, but large proportion of TP53 mutations. Many copy number aberrations and upregulation of interferon response. pMMR. High immune population (B cells and T cells in particular)\nIntClust 5(n=34 tumours):\nMostly SA tumours with high mutational burden. Few copy number aberrations and many up-regulated pathways though. High B cells proportion.\nIntClust 6 (n=20 tumours):\ndMMR tumours with high mutational burden. Upregulation of MYC pathway.\nIntClust 7 (n=19 tumours):\nSA tumours mostly with lower mutational burden, however representation of mutations in NOTCH1 and HRAS. Copy number aberrations are absent. Down-regulation of interferon pathways. High proportion of endothelial cells and fibroblasts.\nIntClust 8 (n=23 tumours):\nMostly SM and SA tumours, low mutational burden but HRAS and ZNF750 mutations. pMMR and many pathways down regulated.\nIntClust 9 (n=10 tumours):\nGroup comprised exclusively of SC-O tumours, pMMR and low tumour burden. However, almost all of them have mutations in TP53 and ZNF750. Upregulation of E2F and G2M checkpoint pathways. Lower immune component. Many copy number aberrations.\nWe also compared this clustering with the iCluster analysis that we presented in our previous version of the manuscript. We saw a moderate agreement between the two partitions (adjusted Rand index=0.43):\n1 2 3 4 5 6 7 8 9\n1 0 11 4 1 10 4 0 0 0\n2 0 0 1 1 2 0 17 0 0\n3 1 0 0 1 0 7 1 7 0\n4 0 0 0 8 0 0 0 3 0\n5 0 0 0 0 0 1 0 13 0\n6 1 0 0 1 0 0 0 0 10\n7 0 0 0 2 22 0 0 0 0\n8 11 2 8 0 0 0 1 0 0\n9 0 0 0 0 0 8 0 0 0\n10 0 2 0 9 0 0 0 0 0\nI should note, even if the authors are able to address the concerns described above, and although I appreciate these exploratory types of integrative analysis, I would argue they cannot replace the more well-established analyses that many readers would expect, if only for clarity and easier interpretation of the results.\nWe agree with the reviewer that this is a complementary analysis that offers an alternate view of the biological features of the disease. In this study we use many different pieces of information that need integration, and we thank the reviewer for the suggestion of comparing this analysis with a purely transcriptomic clustering approach.\nIn the revised manuscript we provide the expression clustering analysis described below.\nIn other words:\n[4.1] How many (purely) transcriptional groups seem to be there? If the exact number is hard to determine, the smallest number that make sense (visually) should be the starting point.\nUsing the most 5,000 variable genes, the Euclidean distance as the dissimilarity matrix and the Ward hierarchical clustering method on the z-scale gene expression matrix, we identified 6 transcriptomic clusters. The optimal number of clusters was selected based on the silhouette measure. We tried several options for the dissimilarity and clustering method and found quite consistent results. The following heatmap shows the expression patterns of the six groups.\nWe provide these data and the associated matrix in the revised manuscript (in Github).\n[4.2] Once these groups have been identified, do they correspond to the ST subtypes (SA vs. SC-E vs. SC-O)?\nThere is not a 1 to 1 correspondence between subtypes and transcriptomic clusters, but there is an association. The following table and the plot, representing the residuals of the Chi-square association test highlights the subtypes more enriched in each group. Exp-Cluster 1 is enriched in OSC (or SC-O) tumours, Exp-Cluster 3 in S tumours and Exp-Cluster 4 in SA tumours. Most of EOSC (or SC-E) tumours are in Exp-Cluster 5.\n1 2 3 4 5 6\nEOSC 3 1 3 3 16 1\nOSC 10 0 0 0 7 0\nS 0 13 15 7 14 14\nSA 0 16 6 32 6 3\n[4.3] Are they enriched for specific mutations or CNAs?\nFor each Exp-Cluster, we run a logistic regression comparing the probability of observing a mutation (from our subset of 84 genes) in the tumours of that Exp-Cluster with the probability of observing that mutation in the rest of the tumours. After multiple test FDR correction, we represented the significant associations in the following plot. Exp-Cluster 1 is characterised by mutations in TP53, ZNF750 and RB1 (although low mutational burden), Exp-Cluster 2 by MSH6 mutations and lack of TP53, Exp-Cluster 3 and 4 by high mutational burden, and Exp-Cluster 6 by lower mutational rates (see figure attached).\nWith respect to copy number aberrations, there are also important differences amongst groups, with group 1 having a very specific copy number pattern, and other groups with very few aberrations. The next plot shows the differences.\n[4.4] Which immune and other non -tumor cell types are present in each of these bulk tumors? How does this correlate with tumor purity?\nUsing EPIC for cell deconvolution, we estimated proportion of different cell types. As the next plot shows, cluster 2 has a high proportion of endothelial cells, cluster 5 (with up-regulation in interferon pathways, see next section) has a high proportion of B cells, CD4+ T cells and macrophages. Cluster 6 has a slightly higher proportion of CD8+ T cells.\n[4.5] Which pathways or genes seems differentially expressed between these groups? Could this be confounded by tumor purity or non-tumor cells or is there evidence that they could represent genuine differences between the cancerous or neoplastic cells?\nWe performed Gene set enrichment analysis on the full expression matrix on the 50 Hallmark Gene sets, and after FDR correction we identified the following enriched pathways in each Exp-Cluster (red means positive normalised enrichment score (NES), and blue negative NES)\nWith respect to tumour cellularity, we had purity estimates obtained from the ASCAT algorithm on the DNA sample for 159 samples. We must note that this method is far from perfect and it produces several 100% purity values in 35 of the cases. However, it should give a rough idea of the association with our clusters. As the plot below shows, in the case of the purely transcriptomic clusters, only Exp-Cluster 5 is composed of lower purity tumours (although with a non-significant adjusted t-test from a linear model, p-value=0.055). In the case of the MOFA IntClusts, we also see a group with lower tumour purity, IntClust 4, but again a non statistically difference in the average purity. Overall, purity does not seem to impact any of the clusters.\nIn summary:\nExp-Cluster 1 is characterised by mainly SC-O tumours and a few SC-E, low mutational burden, but selected mutations in TP53, RB1 and ZNF750. Proficient MMR and pathway activation of E2F targets, G2M checkpoint and MYC targets.\nExp-Cluster 2 is characterised by SM and SA subtypes, 70% of the tumours are MMR deficient and 33% and have several overexpressed pathways, such as Bile Acid Metabolism, myogenesis and UV response down.\nExp-Cluster 3 is enriched in SM tumours, with high mutational burden and almost all of them (96%) MMR deficient (several with MLH1 mutations). They have up-regulation of the MYC pathway and down-regulation of many other pathways.\nExp-Cluster 4 is mostly formed by SA tumours, carrying a large number of mutations (81% of them MMR deficient). Up-regulation of many pathways, such as Estrogen response, inflammatory response, MTORC1 signalling and P53.\nExp-Cluster 5 is mostly characterised by SC-E tumours and activation in the interferon gamma and alpha response pathways. It has a high proportion of B cells, CD4+ T cells and macrophages\nExp-Cluster 6 is mainly formed by SM tumours, with low mutational burden and MMR proficient (94%), with down-regulation of most pathways.\nComparing the purely transcriptomic clusters and the integrative clusters, we see some differences, but the integrative clustering is able to identify interesting tumours with unique combination of mutational and expression profiles, such as IntClust 4.\nClinical data:\n[5] Supplementary Table 1 is a somewhat ambiguous and includes a mix of patient-related and sample-related data with several issues. I recommend that the authors provide two separate tables: one for the patient-related data and one for the sample-related data, each with appropriate IDs (patient and sample).\nWe thank Reviewer2 for sharing his concern with us about Supplementary Table 1. We decided to keep the original and add two supplementary sheets to it following your suggestions, so the reader can select the one that is more appropriate (please see Supplementary Table 1).\n[5.1] Some fields lack sufficient description in the legend, such as \"Age,\" which does not specify whether it refers to the age at diagnosis or at sample procurement. Similarly, “Recurrence”, “TimeAfterInitialDiagnosis_Rec(months)”, “Metastasis” and “TimeAfterInitialDiagnosis_Met(months)” also needs some explanation.\nWe thank Reviewer2 for pointing this out. We corrected this (please see Supplementary Table 1).\n[5.2] Column AE has no name. It contains “Father, Brother” for sample “PD42477a”, but “FamilyCancerLynch” says “no” for this sample.\nWe thank Reviewer2 for pointing out these mistakes, which we fixed (please see Supplementary Table 1).\n[5.3] Some field names differ between the table and the legend. For example, \"MMRstatus\" is mentioned in the legend, but \"MMR.d.MMR.p_simplified\" appears in the table. Same for some IHC fields.\nWe thank Reviewer2 for highlighting these mistakes, which we fixed.\n[5.4] The MMR status is not written correctly for PD42171a (MMR-p instead of pMMR).\nWe thank Reviewer2 for pointing out this mistake, which we have fixed.\n[5.5] For \"MLH1Hypermathylation,\" a mix of blank and \"not tested\" is used. Please use a unique value (either blank or \"not tested\") for all cases to maintain consistency.\nWe thank Reviewer2 and have amended the table as suggested.\n[5.6] Typos: “Hypermathylation”, “Supplemntary”\nWe thank Reviewer2 for pointing out these mistakes, which we have fixed.\n[5.7] There are a number of issues in this table. It would be ideal to ensure that all fields used in the analyses are clearly documented and are error-free. I would recommend fields that cannot be used in any reasonable analyses, because of incomplete or ambiguous data, should be dropped.\nWe thank Reviewer2 for sharing their concerns with us. We have replaced whitespace with NA. Since we have collected all the data variables included, we feel that we are future proofing the data set by including it.\nGermline alleles in sebaceous tumour patients:\n[6.1] This section could be further developed. Associated Figure 6 contains 27 genes, but only a handful of these genes are mentioned in the text. I would recommend improving figure (e.g., removing genes not discussed and adding other plots showing how mutations correlate with phenotypes) or moving it to a Supplementary Figure.\nWe thank Reviewer2 for this recommendation. We have made a new Figure 6 (which is now Figure 7) including only samples mutated for MMR and/or MUTYH genes, as well as clinical details (please see Figure R20/Figure 7). We would like to highlight to the reviewers an important finding. Visual inspection of all germline variants suggested aberrant calls at the POT1 locus which we have removed from the revised figure (below). Indeed, variants found in POT1 were removed as after manual inspection they appeared to be artefactual (please see ‘Germline variant calling’ in the Methods).\nFigure R20/Figure 7 - Germline variant landscape of sebaceous tumour patients. Oncoplot showing the variants present in the normal samples. Constitutional (germline) variants or those with disruptive consequences (i.e. consequences classified as frameshift, nonsense or splice site by VEP, or variants defined as ‘pathogenic’/’likely-pathogenic’ by ClinVar) are shown. Only genes included in the NHS cancer national genomic list for testing and mutated in at least two 26 patients are shown. Abbreviations: P, variant defined as ‘pathogenic’/’likely-pathogenic’ by ClinVar; S, variant identified as significantly depleted by SGE.\n[6.2] VAF and other important information is missing from Supp. Table 12. While the method section states that germline variants were identified with GATK HaplotypeCaller and GenotypeGVCFs, Supp. Table 12 includes fields such as “Tumor_Seq_Allele2” which are typically associated with somatic variant calling. Could the authors clarify this?\nWe thank Reviewer2 for reporting this issue. We have added to the Supp Table 12 the variant allele frequency for the germline analysis (VAF_norm), as well as modified the name of the columns (Tumor_Seq_Allele2 became ALT and Reference_Allele became REF).\nCopy number landscape of sebaceous tumours:\n[7.1] This section could be better developed. It would be ideal if both Figures 4 and 5 could be simplified to include only the key panels that are clearly highlighted to emphasize points made in the main text.\nWe thank Reviewer2 for this comment. We have added titles to make the figures more comprehensive for the reader (please see Figure R21/Figure 5).\nFigure R21/Figure 5 - Somatic copy number alterations in sebaceous tumours. Penetrance plots of somatic copy number alterations in (A) sebaceous adenomas (SA), (B) sebaceoma (SM), (C) extra-ocular sebaceous carcinoma (SC-E), (D) ocular sebaceous carcinoma (SC-O), (E) mismatch repair proficient (pMMR) sebaceous tumours (STs), (F), mismatch repair deficient (dMMR) STs showing loss of MSH2/MSH6, and (G) dMMR STs showing loss of MLH1/PMS2. Frequency plots of copy neutral loss of heterozygosity (cnLOH) in (H) pMMR STs, (I) dMMR STs showing loss of MSH2/MSH6, and (J) dMMR STs showing loss of MLH1/PMS2. Abbreviations: dMMR-MSH2/MSH6, deficient in MMR with MSH2/MSH6 loss of expression by immunohistochemistry; dMMR-MLH1/PMS2, deficient in MMR with MLH1/PMS2 loss of expression by immunohistochemistry; Mb, megabase; min, minimum.\nFigure R1/Figure 6 - Significant focal somatic copy number alterations in sebaceous tumours. Representation of the GISTIC2.0-identified by ST subtype significant focal peak regions (≤ half a chromosome arm) in (A) SA, (B) SM, (C) SC-E) and (D) SC-O. Chromosomal location of the peaks are indicated by the cytoband, with amplifications (gains) shown in red and deletions (losses) shown in blue; cytobands in grey are peaks defined as non-significant. The presence of known tumour suppressor genes or oncogenes within these peak regions are labelled. For the G-score, positive values refer to gains, while negative values are related to deletions. (E) Oncoplot showing the individual tumour samples with the amplification and deletion regions shown in A-D. Samples with no significant gains or losses are indicated by ‘NS’. Five samples were not included in the GISTIC2.0 analysis (see Methods), so no data is shown (indicated by ‘No’). Note: *focal peaks were considered as significant when reaching the 4 criteria: (i) GISTIC2 residual q-value < 0.1, (ii) peak size >= 100,000bp, (iii) concordance with ASCAT CN calls ≥ 0.75, (iv) portion of peak not overlapping with Genome In A Bottle ‘difficult regions’ by >40% (Supplementary Table 9); **approximate cytobands are shown with related genomic coordinates available in Supplementary Table 9.\nNeoantigen and ICI response prediction:\n[8] The last section of the paper (lines 413-482) reports results based on a number of software tools (CODEFACS, TIMER2.0, EPIC, LIRICS, CIBERSORT, pVACseq, and others) and several external datasets that I evaluate as highly speculative. There are insufficient details provided to properly understand the methodology and the results, which are supported by a single, not particularly informative, main figure. At this stage, I would recommend removing this section to ensure clarity and rigor of the whole exome and transcriptome analysis of the paper, or significantly revise this section with more clear details on how the analysis was performed. However, I should note I anticipate the results will likely still be speculative and omission of the ligand-receptor analysis should be considered if validation by IHC/IF staining studies are not being included.\nWe thank Reviewer2 for sharing their concerns about the neoantigen and ligand-receptor analyses. We agree these analyses are purely speculative and require follow-up with further functional analysis. We have largely removed these analyses or moved the remaining elements to the supplementary information with a note that they require validation.\nOther:\n[9.1] There is some inconsistent naming and abbreviation usage for the tumor subtypes across the main text, figures, and their legends. For instance, \"SM\" is used for sebaceoma in Figure 1 and Table S1, while only \"S\" is used in Figure 3 and S5. Similarly, \"PeriOcular,\" \"Peri-ocular,\" and \"ocular\" are all used to refer to ocular sebaceous carcinoma. Please standardize the terminology and abbreviations for each disease subtype.\nWe thank Reviewer2 for pointing out these mistakes, which we have fixed.\n[9.2] According to ST1, sample PD43886c and PD48829c have been used in DNA analysis (DNA_Analysis == Yes). However, I couldn’t find any mutations associated with these samples in ST4.\nWe thank Reviewer2 for pointing out this discordance. We have fixed the table to make it clearer.\n[9.3] SF1A has no legend describing the meaning of each color and does not seem appropriate to clearly support the claims.\nWe thank Reviewer2 for pointing out and we have corrected the figure legend accordingly.\n[9.4] The cohort size is relatively large, considering the rarity of these tumors. However the SC sub-group is quite small, particularly the SC-O subgroup which may limit the power of the analyses. More discussion of this (and other limitations) would be appreciated in a short, dedicated section.\nWe agree that the SC-O collection is small but it remains the biggest collection of this rare tumour type ever studied. In future, additional SC-O may be sequenced and combined with our dataset. We inform the reader about the challenges posed by the SC-O sample number in the results section (“RESULTS - Sample ascertainment and clinico-histopathological features” please see revised Manuscript).\n[9.5] Suggestion: Is it necessary to note that no recurrence or metastasis arose in SA or SM (Line.149) as these are benign tumors that do not usually (or at all?) metastasize. If they do, this should be noted in the text.\nWe thank Reviewer2 and removed the sentence.\n[9.6] Suggestion: It would be ideal to list the 6 countries that these tumors come from (line 81) and mention in discussion any potential biases from tumor selection based on ethnicity or geographical location – for example, authors mention that SC-O is more common in people of Asian descent (40-60% of eyelid tumors). Have the authors investigated the ethnic makeup of their cohort? Is it representative of the global or regional prevalence/distribution of these tumor subtypes? If not, how might this affect the results of their study?\nWe thank Reviewer2 for making this point and have specified the different countries where the tumours came from. An ancestry analysis was performed to explore the genetic origin of the patients. We projected their genotypes against the 1000 genomes dataset (Supplementary Figure 1) revealing largely European and Asian descent as expected.\nReviewer #3 (Remarks to the Author): expertise in neoantigen prediction\nAuthors investigated the genetics and mutational signatures of sebaceous tumors and reached into the conclusion that there are multiple subtypes present. The study also has an interesting observation that some of these tumors also have lynch signatures. The lynch signatures should be investigated further using current non tumor and historical patient samples to reach a definitive conclusion that these mutations were not acquired later in development.\nWe thank Reviewer3 for their comment however, it remains unclear for us what they want us to do. We have performed a germline analysis where we identified known disruptive variants across established genes from the Cancer national germline genomic testing panel performed by the NHS/UK. We observed that several patients carried an inherited mutations in MMR genes and MUTYH as previously described making the aetiology of the Lynch signatures in these patients most likely related to these germline alleles.\nFigure R20/Figure 7 - Germline variant landscape of sebaceous tumour patients. Oncoplot showing the variants present in the normal samples. Constitutional (germline) variants or those with disruptive consequences (i.e. consequences classified as frameshift, nonsense or splice site by VEP, or variants defined as ‘pathogenic’/’likely-pathogenic’ by ClinVar) are shown. Only genes included in the NHS cancer national genomic list for testing and mutated in at least two patients are shown. Abbreviations: P, variant defined as ‘pathogenic’/’likely-pathogenic’ by ClinVar.\nUtilizing some published material to support Sebaceous tumor prevalence in Lynch syndrome would also be informative.\nSebaceous tumours develop in less than 5% of patients with Lynch syndrome (LS) or hereditary non-polyposis colorectal cancer (HNPCC) syndrome, which has also been referred to as Muir-Torre syndrome (MTS) (Marcoval J et al. 2019, Schwartz RA et al. 1995). The incidence of LS/MTS in patients with sebaceous neoplasms ranges between 14 and 50% (Marcoval J et al. 2019, Ponti G et al. 2025, Roberts ME et al. 2013, Roberts ME et al. 2014, Lee BA et al. 2012, Everett JN et al. 2014) suggesting patient ascertainment is an important factor in reporting the frequency of the disease.\nAuthors utilize available pipelines to perform a lot of their computational analysis, it would be helpful if same pipelines are used on some examples of published historical skin tumors such as Melanoma, squamous carcinoma, UV mediated Merkel cell carcinoma for comparison in these figures such that the reader could make comparisons of mutational burden etc. these results may be skewed by the pipelines used.\nWe thank Reviewer3 for they comment. We have extracted data from MSK-impact via the cBioPortal and compare the mutation frequency to our ST data as shown in the plot below.\nFigure R22 - Comparison of mutation count in target regions from four MSK-Impact skin tumour types (Merkel cell [MCC], melanoma, BCC and SCC) to the sebaceous tumours (STs) in this study. The ST cohort are shown all together (red dots), or separated by subtype and mismatch repair (MMR) status (dMMR as green dotsfont and pMMR as black dots). Abbreviations: BCC, basal cell carcinoma; MCC, Merkel cell carcinoma; SCC, squamous cell carcinoma.\nComparing the same panel of mutated genes, we can observe dMMR STs tumours (benign and malignant STs; highlighted in green) have a higher mutation burden compared to other malignant skin tumour subtypes.\nAuthors make some suggestions about the immunogenicity of RCC tumors, this is still currently a puzzle since we don't know what makes RCC tumors highly immunogenic. First the Mutational burden for sure does not explain RCC immunogenicity and secondly the CD8 T cell infiltration in RCC tumors is actually deleterious. So the discussion and comparison to these tumors would need to be carefully worded and considered in the text. Authors make conclusions about the immunogenicity of the tumors. This part requires additional work. First the authors again utilize some pipeline to do their analysis, this would be valid if there are some controls presented for comparison. So the Neoantigen score for Melanomas, squamous carcinomas and UV induced merkel cell carcinoma utilizing published data should be presented as a control comparison.\nWe thank Reviewer3 for their comment about the neoantigen and deconvolution analyses. We have deleted most of these analyses and/or moved them to the supplementary material. We have removed the direct comparisons across tumour types. We feel we present a meaningful analysis of these data.\nAdditionally the authors present histology data in figure 1. As a supplement an immune score calculation utilizing CD4/CD8 infiltration should be presented for each tumor since the authors have access to FFPE blocks. This would make the interpretation of the all neoantigen loads the authors are presenting and would have to validate some of the observations\nWe would not know how to interpret the data based on this as, for example, Brisk and non-Brisk infiltration is more important in melanoma. Again, we moved all the Figures related to these analyses to supplementary material.\nWe thank Reviewer 3 for their very thoughtful comments."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors did a thorough job in addressing the reviewers' concerns and the manuscript has been significantly improved. My comments have been adequately addressed.", "reviewer_2": "The authors’ have made considerable efforts to revise their manuscript and respond to suggested revisions and comments. However, two key issues have not been appropriately addressed and require further attention.\n1. Driver Gene Analysis & False Positives\nIn the initial review, concerns were raised about the likely high proportion of false positives due to the authors' reliance on a single tool (dNdScv) without proper filtering or orthogonal validation. Specifically, SPANXB1 was highlighted as an example of a likely false positive due to the low VAF of its mutations. In response, the authors argue that, based on the interquartile range method, SPANXB1’s VAF is not an outlier. Unfortunately, the authors response misses the point of the major concern raised—that further examination of the significantly mutated gene list is needed to identify and remove/flag likely false positives before publication. For example, in re-analyzing the data provided by the authors, 85% of SPANXB1 mutations have a VAF < 0.25, which is directly at odds with their own Figure R8, where they claim LOH is observed in nearly 70% of cases. Additionally, it is important to note that SPANXB1 is located on the X chromosome, and most of its reported mutations occur in males—where there is no heterozygosity to lose. This raises concerns regarding the LOH claim in Figure R8 and the robustness of the significant mutated gene list.\nImportantly, the suggested additional analyses —VAF distribution, mutation signatures, gene-level RNA-Seq validation (e.g. RNA-VAF vs DNA-VAF), and co-occurrence of LOH with loss-of-function mutations —were essential recommendations for validation/filtering steps to refine the driver mutation gene list before publication (i.e. identify and remove/flag the likely false positives). These analyses should be included in the paper to demonstrate that the final significantly mutated gene list is robust and biologically meaningful.\nOne last additional point, when showing the distribution of mutation consequences by gene (Figure R9), synonymous mutations must be included, as they are essential for interpreting selection patterns (e.g. a high proportion of synonymous mutations suggests weak or no selection).\n2. Integrative Clustering Analysis Based on Unreliable Data\nAnother major concern was that the integrative clustering analysis relies on mutation data dominated by passenger mutations (not just a few as claimed by the authors). The authors did perform a more rigorous expression-based analysis for in their revised manuscript. However, instead of highlighting this new and improved analysis, they relegated it to the supplementary material, while keeping the original problematic figure in the main text. I suggest they do the opposite: put their new analysis in the main text and relegate the one based on passenger mutations in the supp mat (or consider discarding it).\nConclusion\nWe encourage the authors to directly address our original concerns and incorporate the necessary analyses to ensure their findings are well-supported.", "reviewer_4": "I believe the authors have appropriately addressed all of the valid concerns raised by Reviewer #3. Regarding the previous question about 'these results may be skewed by the pipelines used' and whether 'it would be helpful if the same pipelines were applied to some examples of published historical skin tumors,' while this is a valid concern regarding potential bias arising from differences in bioinformatics pipelines, however, re-analyzing the data using the same pipeline would not eliminate systematic biases between studies. These biases also arise from variations in sequencing methods and sample preparation protocols used by different research institutes. Additionally, publicly available datasets, such as MSK-IMPACT, utilize custom panels with different targeted regions, re-analyzing these datasets using the same bioinformatics pipeline would not resolve the inherent systematic differences."}, "author_response": "> We thank Reviewer #1 for reviewing this paper and their insights into our study.\nReviewer #2/3 (Remarks to the Author):\nThe authors’ have made considerable efforts to revise their manuscript and respond to suggested revisions and comments. However, two key issues have not been appropriately addressed and require further attention.\n> We thank Reviewer #2 for reviewing our revised manuscript and apologise if we misunderstood some of their points.\n1. Driver Gene Analysis & False Positives\nIn the initial review, concerns were raised about the likely high proportion of false positives due to the authors' reliance on a single tool (dNdScv) without proper filtering or orthogonal validation. Specifically, SPANXB1 was highlighted as an example of a likely false positive due to the low VAF of its mutations. In response, the authors argue that, based on the interquartile range method, SPANXB1’s VAF is not an outlier. Unfortunately, the authors response misses the point of the major concern raised—that further examination of the significantly mutated gene list is needed to identify and remove/flag likely false positives before publication. For example, in re-analyzing the data provided by the authors, 85% of SPANXB1 mutations have a VAF < 0.25, which is directly at odds with their own Figure R8, where they claim LOH is observed in nearly 70% of cases. Additionally, it is important to note that SPANXB1 is located on the X chromosome, and most of its reported mutations occur in males—where there is no heterozygosity to lose. This raises concerns regarding the LOH claim in Figure R8 and the robustness of the significant mutated gene list.\n> We thank Reviewer #2 for noting that SPANXB1 is the only driver genes located on Chr X which could cause misinterpretation in the copy number analysis and more specifically in cnLOH. In agreement with Reviewer #2, we then decided to remove SPANXB1 from Figure R8 which is included in Supplementary Figure 4.\nSupplementary Figure 4 – Driver genes and associated changes. (A) Table showing the proportion of somatic mutations also present in RNAseq. (B) Barplot showing the proportion of the different variant types across driver genes. (C) Barplot showing the proportion of the somatic copy number alterations (SCNAs) across autosomal driver genes. (D) Scatter plot showing missense versus truncating mutation profiles across the driver genes. (E) Barplot showing the distribution of the variant allele frequency (VAF) across the driver genes. Each dot is a tumour. Of note, SPANXB1 is located on the X chromosome which is likely reflective of the low VAF we observe. (E) Tileplot showing the alteration of the driver genes across sebaceous tumours and pan-cancer analyses. Genes are from Bailey MH et al., 17. The gene role is taken as defined by Bailey MH et al., in 17 and COSMIC Cancer Gene Census database for TCGA cancers and sebaceous tumours, respectively. Abbreviations: NA, not defined; onco, oncogene; onco.tsg, oncogene and tumour suppressor genes; tsg, tumour suppressor gene.\nImportantly, the suggested additional analyses —VAF distribution, mutation signatures, gene-level RNA-Seq validation (e.g. RNA-VAF vs DNA-VAF), and co-occurrence of LOH with loss-of-function mutations —were essential recommendations for validation/filtering steps to refine the driver mutation gene list before publication (i.e. identify and remove/flag the likely false positives). These analyses should be included in the paper to demonstrate that the final significantly mutated gene list is robust and biologically meaningful.\n> We thank Reviewer #2/3 for this suggestion and we have added Figure R5 and the table comparing the number of mutations found in the DNA and RNA data, which reinforce the fact that the data are of high quality, to Supplementary Figure 4.\nSupplementary Figure 4 – Driver genes and associated changes. (A) Table showing the proportion of somatic mutations also present in RNAseq. (B) Barplot showing the proportion of the different variant types across driver genes. (C) Barplot showing the proportion of the somatic copy number alterations (SCNAs) across autosomal driver genes. (D) Scatter plot showing missense versus truncating mutation profiles across the driver genes. (E) Barplot showing the distribution of the variant allele frequency (VAF) across the driver genes. Each dot is a tumour. Of note, SPANXB1 is located on the X chromosome which is likely reflective of the low VAF we observe. (E) Tileplot showing the alteration of the driver genes across sebaceous tumours and pan-cancer analyses. Genes are from Bailey MH et al., 17. The gene role is taken as defined by Bailey MH et al., in 17 and COSMIC Cancer Gene Census database for TCGA cancers and sebaceous tumours, respectively. Abbreviations: NA, not defined; onco, oncogene; onco.tsg, oncogene and tumour suppressor genes; tsg, tumour suppressor gene.\nOne last additional point, when showing the distribution of mutation consequences by gene (Figure R9), synonymous mutations must be included, as they are essential for interpreting selection patterns (e.g. a high proportion of synonymous mutations suggests weak or no selection).\n> We thank Reviewer #2/3 for this suggestion and we have modified Figure R9 by including synonymous (silent) variants which we include into Supplementary Figure 4.\nAs shown, synonymous mutations are not major events in these driver genes.\nSupplementary Figure 4 – Driver genes and associated changes. (A) Table showing the proportion of somatic mutations also present in RNAseq. (B) Barplot showing the proportion of the different variant types across driver genes. (C) Barplot showing the proportion of the somatic copy number alterations (SCNAs) across autosomal driver genes. (D) Scatter plot showing missense versus truncating mutation profiles across the driver genes. (E) Barplot showing the distribution of the variant allele frequency (VAF) across the driver genes. Each dot is a tumour. Of note, SPANXB1 is located on the X chromosome which is likely reflective of the low VAF we observe. (E) Tileplot showing the alteration of the driver genes across sebaceous tumours and pan-cancer analyses. Genes are from Bailey MH et al., 17. The gene role is taken as defined by Bailey MH et al., in 17 and COSMIC Cancer Gene Census database for TCGA cancers and sebaceous tumours, respectively. Abbreviations: NA, not defined; onco, oncogene; onco.tsg, oncogene and tumour suppressor genes; tsg, tumour suppressor gene.\n2. Integrative Clustering Analysis Based on Unreliable Data\nAnother major concern was that the integrative clustering analysis relies on mutation data dominated by passenger mutations (not just a few as claimed by the authors). The authors did perform a more rigorous expression-based analysis for in their revised manuscript. However, instead of highlighting this new and improved analysis, they relegated it to the supplementary material, while keeping the original problematic figure in the main text. I suggest they do the opposite: put their new analysis in the main text and relegate the one based on passenger mutations in the supp mat (or consider discarding it).\n> We disagree with this comment as we think there is a misunderstanding in the goals and the details of our Integrative Clustering. Methods such as MOFA or IntClust+, the ones we performed on this manuscript, do not simply perform clustering of tumours using the features as inputs, as the hierarchical clustering based on expression data mentioned by the reviewer would do. Instead, they fit a version of a factor analysis model that includes regularisation on the parameters. This means that the model “shuts off” mutations or expression values that occur randomly in tumours and not systematically in subsets of tumours.\nIn fact, both MOFA and IntClust were designed as methods to identify driver genes, and have been used for that task in many studies, such as identification of driver genes in triple negative breast cancer tumours (https://doi.org/10.1016/j.celrep.2023.113564), identification of key biological processes in kidney cancer (doi: 10.1038/s44320-024-00072-3), discovery of driver-based breast cancer subtypes (doi:10.1038/nature10983), driver genes discovery in lung cancer (doi: 10.1038/s41467-019-11276-9), or Crohn disease (doi: 10.1093/ibd/izaa281).\nOur integrative analysis identified several driver processes, one of them characterised by overall mutation burden, and others by only a handful of mutations, typically ZNF750, TP53, RREB1 and NOTCH2. These are the genes with high values of the parameters that will contribute to the final Integrative Clustering of tumours. It is true that there can be passenger genes within the set of input features, but MOFA reduces (or sets to zero) the coefficient of these genes, having no impact in the integrative clustering. We think that our analysis shows the interplay between mutations and gene expression, providing a richer view of the driver event landscape than a simple clustering based on gene expression. Comparing both approaches, as we showed in our previous response, our integrative clustering provides additional insights.\nIn summary, the integrative clustering is not dominated by passenger mutations. MOFA identifies a handful of driver genes and uses only those to perform the clustering of tumours. \nConclusion\nWe encourage the authors to directly address our original concerns and incorporate the necessary analyses to ensure their findings are well-supported.\nReviewer #4 (Remarks to the Author):\nI believe the authors have appropriately addressed all of the valid concerns raised by Reviewer #3. Regarding the previous question about 'these results may be skewed by the pipelines used' and whether 'it would be helpful if the same pipelines were applied to some examples of published historical skin tumors,' while this is a valid concern regarding potential bias arising from differences in bioinformatics pipelines, however, re-analyzing the data using the same pipeline would not eliminate systematic biases between studies. These biases also arise from variations in sequencing methods and sample preparation protocols used by different research institutes. Additionally, publicly available datasets, such as MSK-IMPACT, utilize custom panels with different targeted regions, re-analyzing these datasets using the same bioinformatics pipeline would not resolve the inherent systematic differences.\nWe thank Reviewer #4 for raising this important point, and we agree. Quantitative comparisons between sample groups must be treated objectively, particularly when comparing cohorts with both high and low mutational load tumors, and when data are generated on different platforms. We have therefore included a caveat in the discussion of these data and discuss them qualitatively."}, {"round": 3, "reviewer_comments": {"reviewer_2": "Despite previous requests for more robust validation of the nominated drivers, the authors persist with their original list of 28 drivers, relying solely on the outputs from dNdScv. While dNdScv is effective, it is prone to reporting false positives especially in high TMB contexts. This is concerning as many of the nominated drivers are not established cancer-associated genes or have not been identified in large pan-cancer driver mutation analysis (for example nearly half (13/28) not recognized in a characterization of cancer driver genes in the TCGA PanCancer analysis by Bailey et al). The authors do mention using oncodriveFML as a complementary tool in their methods section, yet these results are never discussed in the manuscript, nor are they integrated with those from dNdScv. The raw dNdScv output is still not included in this revision as originally requested\nDespite clear evidence that indicates SPANXB1 is a false positive, the authors continue to nominate it as a driver gene. They incorrectly attribute the unusually low variant allele frequency (VAF) of its mutations to its location on the X chromosome in their revised manuscript. Indeed, VAF of mutations on the X chromosome should be twice that of heterozygous mutations on autosomes in males.\nIn Supplementary Table 5B, where the results of a driver analysis were supposedly restricted to drivers identified by Bailey et al., genes such as TTN, OBSCN, and NEB with extremely low p-values were included. These genes are known for accumulating passenger mutations due to their length and are not listed among the drivers in Bailey et al. study. This discrepancy indicates that either the analysis was not conducted properly or the incorrect results were presented.\nFig S4A contains obvious errors. How could 676 SNV mutations be observed with >10X coverage if the total number of SNV mutations is 667? The authors also missed the point since the idea here was to show that most mutations in *each* of the nominated drivers are observed in the RNA-seq, in contrast to passenger mutations in other genes that are often not expressed. As is, the table does not inform about the individual plausibility of the nominated drivers. Figure S4C and S4E also miss the point and do not inform the plausibility of the nominated drivers. The idea was to show if mutations often co-occur with cnLOH or cnLoss in nominated TSGs (and hence have higher VAF), in contrast to passenger mutations in *other* genes, as would be expected. Additionally, in Figures S4B and S4C, the y-axis is inaccurately labeled as 'count'.\nThe purpose of these additional analyses we requested was to identify which candidate drivers are substantiated by orthogonal evidence and to pinpoint likely false positives. Unfortunately, the authors chose to present inconsequential data and assert that their initial list of drivers remains valid without requiring any modifications.\nI do still think this is an important resource paper to the cancer genomics community.", "reviewer_5": "Dear Dr McGinnis.\nI have gone through the questions from Reviewer #2 and the authors' responses. Although some of the points are addressed in detail, there are still some points that require clarification.\na) In the author's response, it was mentioned that 15 out of 28 predicted drivers are present in the TCGA PanCancer analysis by Bailey et al., and the remaining ones are linked to tumorigenesis based on literature evidence. I think this is informative, and the authors may consider providing this information in the supplementary table and also mentioning this overlap number in the results section.\nb) The results of restricted hypothesis testing using Bailey et al., (299 cancer driver genes) reported in Supplementary Table 5B still show genes such as TTN, OBSCN, and NEB, which are likely false positives. The authors stated that they were removed, but they're not reflected in the table yet.\nc) The results of oncodriveFML were not directly included in the supplementary table, so it is difficult to compare the results. The GitHub link to the oncodriveFML output is also not working. It would be helpful for the readers to know at least what fraction of the 28 potential drivers predicted by the dNdScv overlap with the oncodriveFML significant hits. In particular, whether the four drivers (NOTCH1, RREB1, HRAS, TP53) used for molecular subtyping are supported by oncodriveFML or not."}, "author_response": "NCOMMS-24-26940B: The molecular cartography of malignant and benign\nsebaceous tumours\nReviewer #2 (Remarks to the Author)\nDespite previous requests for more robust validation of the nominated drivers, the authors persist with their original list of 28 drivers, relying solely on the outputs from dNdScv. While dNdScv is effective, it is prone to reporting false positives especially in high TMB contexts.\nThe original dNdScv publication (Martincorena et al., PMID: 29056346) has been cited over 1,500 times and is widely used in the field, including by the Hartwig Medical Foundation (PMIDs: 29056346, 31645765), the international PCAWG project (PMID: 33831375), and in studies of positive selection in normal tissues, such as during clonal hematopoiesis (PMID: 38744975).\nIndependent benchmarking across multiple cohorts—including exome sequencing datasets—has shown that dNdScv achieves a low false discovery rate and high sensitivity (e.g., PMIDs: 38808071, 31363082). We acknowledge Reviewer 2’s previous comment regarding the potential impact of mutational load on false positives. To address this, we filtered our dataset by mutation count and repeated the driver gene analysis and supplied this in the previous revision of the paper. The results were consistent, confirming that our predictions are robust. While no driver gene method is perfect for samples with high mutational burden, dNdScv performs better than most alternatives due to its use of a background mutation rate–adjusted statistical model.\nThis is concerning as many of the nominated drivers are not established cancer-associated genes or have not been identified in large pan-cancer driver mutation analysis (for example nearly half (13/28) not recognized in a characterization of cancer driver genes in the TCGA PanCancer analysis by Bailey et al).\nReviewer 2 suggests that our list of driver genes includes false positives. We would like to clarify that the Bailey et al. dataset is not a definitive reference or ground truth. The tumours included in that study did not include sebaceous tumours—or, more broadly, adnexal skin tumors. As such, the absence of certain genes from their list may not diminish or invalidate our findings. Notably, the Bailey list also omits several well-established cancer genes, including ETV6, POT1, and ALPK1, the latter of which is a known driver in an adnexal tumour type closely related to sebaceous tumours. Our study was explicitly designed to identify novel driver genes, including those that may not have been previously reported. Filtering or disregarding genes because they do not appear in the Bailey et al. list may risk missing biologically meaningful discoveries. In response to the specific comment about the 28 genes we report, we have annotated the list below. Importantly, 15 of these 28 genes are included in the Bailey et al. list. The other genes are tightly linked to tumorigenesis. We believe this supports the relevance of our findings and underscores the importance of considering tumor-type–specific contexts when evaluating driver gene predictions.\nThe authors do mention using oncodriveFML as a complementary tool in their methods section, yet these results are never discussed in the manuscript, nor are they integrated with those from dNdScv. The raw dNdScv output is still not included in this revision as originally requested.\nWe thank Reviewer 2 for this comment. Firstly, we provide the oncodriveFML analysis in the paper and point readers at this list in the manuscript. We do not discuss it extensively because driver gene analysis is just one part of the paper but, as requested, we provide this analysis as a complement to dNdScv. We do not feel that it is valuable to integrate driver predictions from both methods because there is no clear rubric as to how to do this and there is no way of knowing if we have over filtered the data or introduced errors. The raw dNdScv output was provided with the revision we returned (see NatureComs_dNdScv_Output.xls).\nDespite clear evidence that indicates SPANXB1 is a false positive, the authors continue to nominate it as a driver gene. They incorrectly attribute the unusually low variant allele frequency (VAF) of its mutations to its location on the X chromosome in their revised manuscript. Indeed, VAF of mutations on the X chromosome should be twice that of heterozygous mutations on autosomes in males.\nThere are many genes that have subclonal drivers, for example low VAF EGFR (PMID: 33869038) and NOTCH1 (PMID: 35430388) mutations are common. We cannot exclude a gene from an analysis because we can’t explain how it is contributing to tumorigenesis. We suspect that the reason why this gene has a significant P value is that it mainly has frameshift alterations. It is only called a driver in extra-ocular SC with a borderline P value, but as above we do not feel that it is appropriate to exclude data from an analysis when we have no technical evidence suggesting that the result is not correct. The statement about VAF is only correct in males. Reviewer 2 will note that although SPANXB1 is called as a driver we have now focused our analysis on drivers on the autosomes, so it does not feature in any of the downstream analyses. We hope that this strikes a balance between transparency and pragmatism given Reviewer 2s comment.\nIn Supplementary Table 5B, where the results of a driver analysis were supposedly restricted to drivers identified by Bailey et al., genes such as TTN, OBSCN, and NEB with extremely low p-values were included. These genes are known for accumulating passenger mutations due to their length and are not listed among the drivers in Bailey et al. study. This discrepancy indicates that either the analysis was not conducted properly or the incorrect results were presented.\nWe acknowledge that we made a mistake with this analysis, and we thank Reviewer 2 for noticing this. We have corrected the table and supply it with the revised manuscript.\nFig S4A contains obvious errors. How could 676 SNV mutations be observed with >10X coverage if the total number of SNV mutations is 667? The authors also missed the point since the idea here was to show that most mutations in *each* of the nominated drivers are observed in the RNA-seq, in contrast to passenger mutations in other genes that are often not expressed. As is, the table does not inform about the individual plausibility of the nominated drivers. Figure S4C and S4E also miss the point and do not inform the plausibility of the nominated drivers. The idea was to show if mutations often co-occur with cnLOH or cnLoss in nominated TSGs (and hence have higher VAF), in contrast to passenger mutations in *other* genes, as would be expected. Additionally, in Figures S4B and S4C, the y-axis is inaccurately labeled as 'count'.\nThere are 676 variants that are in driver genes with RNA-Seq >10x. Of these 676 variants, 667 variants (called from the exome data) were re-called (98.7%) in the RNA-Seq data. Given that >98.7% of all SNV mutations validated (as above) we feel that we answered the reviewer’s question as written. i.e. the analysis reveals that we can see the driver events in the deep RNA-Seq data we generated.\nSimilarly, with Figure S4C and E we delivered the analysis requested by the reviewer, so we are unsure of what further to provide.\nPlease note, we are absolutely committed to responding to reviewer 2s comments, but we cannot interpret what they are requesting unless they provide a full description.\nThe purpose of these additional analyses we requested was to identify which candidate drivers are substantiated by orthogonal evidence and to pinpoint likely false positives. Unfortunately, the authors chose to present inconsequential data and assert that their initial list of drivers remains valid without requiring any modifications.\nWe thank you for your attention to detail. We certainly didn’t intend to present data that Reviewer 2 didn’t feel answered the question. We genuinely don’t know what we haven’t done here.\nI do still think this is an important resource paper to the cancer genomics community.\nWe thank Reviewer 2 for their help in improving our manuscript.\nReviewer #5 (Remarks to the Author):\nI have gone through the questions from Reviewer #2 and the authors' responses.\nAlthough some of the points are addressed in detail, there are still some points that\nrequire clarification.\na) In the author's response, it was mentioned that 15 out of 28 predicted drivers are present in the TCGA PanCancer analysis by Bailey et al., and the remaining ones are linked to tumorigenesis based on literature evidence. I think this is informative, and the authors may consider providing this information in the supplementary table and also mentioning this overlap number in the results section.\n> We thank Reviewer #5 for its comment and have amended the manuscript accordingly. The data is in Supplementary Table 5B.\nb) The results of restricted hypothesis testing using Bailey et al., (299 cancer driver genes) reported in Supplementary Table 5B still show genes such as TTN, OBSCN, and NEB, which are likely false positives. The authors stated that they were removed, but they're not reflected in the table yet.\n> These genes are not in the table – after the comments from Reviewer 2 they were removed.\nc) The results of oncodriveFML were not directly included in the supplementary table, so it is difficult to compare the results. The GitHub link to the oncodriveFML output is also not working. It would be helpful for the readers to know at least what fraction of the 28 potential drivers predicted by the dNdScv overlap with the oncodriveFML significant hits. In particular, whether the four drivers (NOTCH1, RREB1, HRAS, TP53) used for molecular subtyping are supported by oncodriveFML or not.\n> We thank Reviewer #5 for this comment. The results of oncodriveFML have been put on the GitHub repository. Out of the 25 driver genes called in the entire cohort by dNdScv, 22 were also significantly found by oncodriveFML (88%) (please see “oncodrivefml_variants-oncodrivefml_dNdScvComparison”). The last 3 genes found in specific groups could not be confirmed as we did not run oncodriveFML that way."}]

Decision letter