Molecular landscapes of oral cancers of unknown etiology

Patient clustering based on de novo mutational signatures

To understand the possible origins and specific molecular features of the NIRF OCs, we performed an in-depth analysis of mutational signatures for non-HPV related HNSCC, ie OC and LX squamous cell carcinoma (SCC).

Mutational signature extraction was conducted on 425 HNSCC samples involving 314 (74%) OC and 111 (26%) LX SCC from The Cancer Genome Atlas (TCGA). After initial sample filtering (see Methods and Supplementary Figure S2), we performed mutational signature analyses of the 347 remaining samples using the 384-channel (i.e. transcription strand-biased) signature versions that ensure a higher sensitivity. We identified four distinct de novo signatures SBS384A-C (Figure 1a) closely matching several known COSMIC signatures, as follows: SBS384A - SBS4 (tobacco smoking exposure); SBS384B - SBS2 and SBS13 (APOBEC enzyme mutagenesis); SBS384C - SBS1 (5-methyl-cytosine deamination), and SBS384D - a set of mutually similar signatures SBS5 (clock-like, unknown etiology), SBS92 (new tobacco smoking signature) and SBS16 (unknown etiology, putatively associated with alcohol consumption) (Supplementary Figure S3).

• Download figure
• Open in new tab

Figure 1:

Mutational signatures that characterize the major OCSCC and LXSCC subsets. (a) De novo transcriptional strand-asymmetry signatures (left panels, signatures SBS384A through SBS384D extracted by NMF, see Methods for details) determine four major patient clusters (heatmap, right panel). The heatmap is based on relative per-sample contribution of each mutational signature. Heatmap annotation tracks and color legends: Mutation count histogram, anatomical subsites (LX = larynx, OC = oral cavity), tobacco smoking status (NS = non-smoker, S = smoker), HPV infection status and no identified risk factors (NIRF). (b) Per-sample and per-cluster distributions of the total mutation counts. (c) Absolute mutation counts attributed to each de novo signature in each cluster. Cluster names reflect the matching putative mutagenic process/COSMIC signature with the highest contribution in each cluster (5mC deamin = spontaneous deamination of methylated cytosine/COSMIC SBS1).

Based on the relative per-sample contribution of the four de novo signatures, the samples were grouped into four distinct patient/signature-based clusters (Figure 1a). In each cluster, the absolute mutation counts were dominated by the activity of one de novo signature (Figure 1b), with more low-mutation burden samples in Clusters 1 (median=139, range=24-338) and 3 (median=230, range=62-457) and high mutation-count samples in Clusters 2 (median=344, range=86-1835) and 4 (median=682, range=153-2193) (Figure 1c). The high relative enrichment of each de novo mutational signature (Figure 1a) reflected well the absolute counts of signature-attributed mutations (Figure 1b). Specifically, signature SBS384C (putative SBS1) appeared to operate at the highest levels in Cluster 1 (renamed 5mC deamin cluster, n=104) (Figure 1b). The mutation counts supported the overall high relative enrichment of SBS384B (putative APOBEC, n=80) in Cluster 2 (renamed APOBEC cluster). SBS384D (resembling COSMIC signatures SBS16, SBS5, or SBS92) was enriched in Cluster 3 (renamed SBS16-like, n=110). Cluster 4 (renamed tobacco smoking cluster, n=53) was dominated by the SBS384A signature (putative SBS4) and by the associated high mutation counts.

Clinico-pathological characteristics match the patient-signature clusters

The NIRF cancer patients were more often found in the 5mC deamin and APOBEC clusters (61% or 32% of 72 NIRF cases, respectively) (Figure 2a). The 5mC deamin cluster further contained the majority of the young patients (< 40 years; 80%, or 8/10), while most of the 91 elderly patients (≥70 years) were equally distributed between the 5mC deamin and APOBEC and clusters (37%, or 34/91) and 35%, or 32/91), respectively)(Figure 2b). The SB16-like and the tobacco smoking clusters included the majority of drinkers (77.8 %, or 14/18 and 91.1 %, or 51/56, respectively) and smokers (100 %, or 53/53 and 95.4 %, or 103/108, respectively) (Figure 2c and 2d). We further observed a higher proportion of female patients in both the APOBEC and 5mC deamin clusters whereas males were more represented in the tobacco-smoking and SBS16-like clusters (Figure 2e).

• Download figure
• Open in new tab

Figure 2:

Mutational signature cluster-based demographic and clinical presentation summaries. (a) Proportion of the NIRF cases across clusters. (b) Distribution of patients by age (young patients ≤40; 40<middle aged<70; old ≥70) per cluster. (c) Relative alcohol consumption and (d) tobacco smoking status distributions across clusters (ND = non-drinker, D = drinker, NS = non-smoker, S = smoker). (e) Relative anatomical sub-site distribution in each cluster, for all analyzed patients. (f) Relative gender-based patient distribution in each cluster. NA status data are not plotted.

The cancer anatomical subsite distribution also tended to follow the clusters (Figure 2f and Supplementary Figure 5). A half (52%, or 49/95) of the cancer originating from the larynx were in the tobacco smoking cluster while tumors affecting the floor of the mouth were enriched in the SBS16-like cluster (63% = 33/54). Moreover, more than half of tumors involving the oral tongue were in the 5mC deamin cluster (54% = 55/103), while cancer subsites involving the hard palate (83 % = 5/6) or the alveolar ridge (75% = 12/16) were found mainly in the APOBEC cluster. Finally, tumors affecting the buccal mucosa were detected in the two NIRF enriched clusters, 5mC deamin (42%) = 8/19) and APOBEC (9 of 19, or 47%).

Collectively, these data show that mutational signature-based clustering recapitulates the clinico-pathological characteristics of patients with or without identified risk factor and supports the potential clinical relevance of the patient stratification according to the de novo mutational signature contents.

COSMIC mutational signatures operating in OC and LX SCC

To further characterize the individual de novo signatures, which can often represent mixed patterns, we performed their further decomposition into the COSMIC reference signatures. This process, performed by using SigProfilerExtractor (SPE), determined 8 COSMIC signatures: SBS1, SBS2, SBS4, SBS5, SBS6, SBS7a/b, SBS13 and SBS16. COSMIC signature SBS4 (tobacco smoking) constituted 60% of the SBS384A de novo signature. Nearly 90% of SBS384B (cosine similarity = 0.989) was jointly generated by COSMIC signatures SBS2/SBS13 (mutagenesis by APOBEC enzymes). SBS384D consisted of a combination of three distinct mutational signatures of approximately equal contribution: COSMIC SBS1, resulting from methylcytosine deamination; SBS6, a mutational process putatively linked to defective DNA mismatch repair, and SBS5, a clock-like signature of unknown origin. A low contribution of SBS7a and SBS7b, linked to UV exposure, was observed in SBS384C (8% and 7%, respectively). Finally, SBS384D was mainly the result of SBS5 (70%), with a contribution of SBS16 (20%), a signature putatively linked to alcohol exposure (Supplemental Figure S6a).

SBS5 can mask mutational signatures associated with exogenous exposures

Sample clustering based on the COSMIC signature decomposition and attribution did not reproduce the sample clustering based on the de novo signature contents (Figure 3a, and Supplemental Figure S6b). This was mainly due to the activity broadly attributed to signature SBS5 (such as its 29.06 % attribution in SBS384C) (Supplemental Figure S6a). The SPE tool assumes that SBS1 and SBS5 are always present and active in all tissues ²⁶. We hypothesized that this could overestimate the respective contributions of SBS1 and SBS5 in the analyzed data set. To test this hypothesis, we performed a signature analysis complementary to the SPE decomposition, by using the Mutational Signature Analysis (MSA) ²⁷ attribution approach (see Methods). Using the signature set proposed by the SPE decomposition (SBS1, SBS2, SBS4, SBS5, SBS6, SBS7a/b, SBS13 and SBS16), we observed an overall reduced attribution of SBS5 resulting in restored stability of patient clusters (Figure 3a, Supplementary Figure S6c). Next, we applied the MSA-based attribution using the entire COSMIC signature set (70 signatures). This approach resulted in cluster stability consistent with the de novo signature-based clustering (Figure 3a, Supplementary Figures S7a,b), and the main contributions of each principal mutational process to a given patient cluster were restored (Figure 3b, Supplementary Figure S7c). Thus, the overfitting of the SBS5 signature within the SPE decomposition step due to the default inclusion of the signature can be effectively resolved by the MSA-based optimized refitting approach.

• Download figure
• Open in new tab

Figure 3:

Attribution of de novo mutational signatures to COSMIC signatures by MSA fitting. (a) Heatmap showing de novo mutational signatures across patients that restore the original clusters described in Fig. 1a. Heatmap tracks show the mutation count histogram, anatomical subsites (LX = larynx, OC = oral cavity), no identified risk factor status (NIRF). Cluster color labels under the heatmap reveal the comparative sample redistribution outcomes after SigProfiler-based fitting and MSA fitting. The MSA-based fitting was performed by using the same signature set identified by SigProfiler fitting, or with a signature set encompassing all COSMIC signatures (v3.2). (b) Bubble plot documenting the predominant, expected or newly revealed signatures in each cluster, with median mutation count attributed to each signature (bubble size) and positive sample frequency (bubble color). (c) Relative distribution of mutational signatures listed in (b) according to the specific head-and-neck anatomical subsites.

The possible masking of other signatures by SBS5 was based on several lines of evidence. Firstly, SBS5 in its transcriptional strand bias version, is highly similar with signatures SBS16 and SBS92 (Supplementary Figure S3). Indeed, MSA analysis could resolve the broad SPE-suggested presence of SBS5 into biologically plausible assignment of a newly described tobacco-smoking-related signature SBS92 almost exclusively to the LX tumors in the SBS4 dominated tobacco-smoking cluster (Figure 3b and Supplementary Figure S7a,c). Secondly, the resolution of SBS5 revealed an increased presence of SBS16 in smokers who were drinkers without compromising the stability of the de novo signature-based clusters (Supplementary Figures S7 and S8).

Tobacco smoking induces varying, site-specific mutational processes

Interestingly, by using the NNLS-based MSA method of signature fitting, signature SBS4 related to tobacco smoking exposure (particularly to the effects of polycyclic aromatic hydrocarbons such as benzo[a]pyrene ^20,28–30) was detected in tumors in the tobacco smoking cluster (involving mostly the LX) but not in the SBS16-like cluster (involving mostly the OC) (Figures 3b, 4 and Supplemental Figure S7a,c) even though both clusters consisted mainly of smokers (Figure 2d). This was further supported by the identification of recently reported tobacco smoking-related signature SBS92 which we observed in the LX group alongside SBS4 (Figure 3b, 4, Supplemental Figure S7a,c). This finding is suggestive of selective, anatomical site-specific mutagenic effects of tobacco smoke, with distinct mutational processes operating in the smokers’ laryngeal tissues (SBS4, SBS92) and in the OC (SBS16).

Tumors in the NIRF clusters are driven by endogenous mutagenic processes

NIRF-enriched clusters were characterized by increased levels of signature SBS1 in terms of both absolute and relative mutation counts (Figure 4 and Supplemental Figure S7a). Moreover, APOBEC signatures SBS2 and SBS13 were predominant in the APOBEC cluster (Figure 4 and Supplemental Figure S7a). Interestingly, in the APOBEC and 5mC deamin clusters, NIRF and drinker/smoker patients showed similar mutation spectra, except for the low presence of transcription strand biased T>C mutation pattern in the smoking/drinking patients’ tumors. The low presence of SBS16 signature in these patients suggests a weak contribution of alcohol/tobacco in the total mutation count in these clusters (Supplemental Figure S9). In sum, NIRF OC SCC are defined by elevated mutagenic activities of endogenous processes that result in the enriched presence of COSMIC signatures SBS1 and SBS2/13.

• Download figure
• Open in new tab

Figure 4:

The main cluster-defining de novo signatures fitted by using MSA to their COSMIC signature contents. Top-left panels: the 96-channel COSMIC signatures; bottom-left panels: the corresponding transcriptional-strand asymmetry signature versions. The boxplots (right panels) show the per-cluster distribution of mutation counts attributed to each COSMIC signature.

Cancer driver gene mutations in the NIRF-enriched clusters

The presence of distinct mutational processes predominantly operating in each cluster suggests potential differences in the natural histories of the disease development. To test this hypothesis, we determined mutated genes that were under selective pressure (see Methods). We identified 26 mostly recurrently mutated genes under positive selection with the following distribution: 5mC deamin cluster – 18/26, APOBEC cluster – 17/26, SBS16-like cluster – 11/26, tobacco smoking cluster – 6/26 (Figure 5a).

• Download figure
• Open in new tab

Figure 5:

Cluster distribution of the identified non-silently mutated genes that are under selective pressure. (a) Summary of the mutated genes under selective pressure identified by the non-synonymous/synonymous mutation ratio approach (see Materials and Methods). The number of samples in each cluster is shown in the top row, and number of cases with a mutated gene (with symbols listed on the left) is shown on the right. P-values for substitutions were obtained by Likelihood-Ratio Tests as described ⁶⁷ and cluster frequency indicate percentage of sample in each cluster mutated for a given gene. (b) Venn diagram showing the distribution of mutated genes under selective pressure across clusters. (c) Relative distribution of signature-matched mutation processes for all mutated genes in each cluster (left) and mutated genes that are under selective pressure (right). Numerical values in the graphs indicate the underlying mutation counts.

The mutated tumor suppressor genes TP53, CDKN2A and FAT1 were selected for in all clusters, whereas various mutated epigenetic modifiers under positive selection exhibited more cluster specific distribution. Mutations in the histone methyltransferase NSD1 were restricted to the tobacco smoking clusters (tobacco smoking and SBS16-like) while KMT2B mutants were only detected in the 5mc deamin cluster. Interestingly, genes involved in mitogen-activated protein kinase cascade components (HRAS and EPHA2) were mutated and selected for only in the NIRF-enriched clusters. The altered antigen presentation system (mutated HLA-A, HLA-B, B2M genes) and immune cytotoxic response (mutated CASP8 gene) were also specific to the NIRF-enriched clusters (APOBEC and 5mC deamin), suggesting the presence of immune pressure in these clusters (Figure 5b).

The functional impact of mutational signatures

To address the impact of the key mutational signatures on all mutated genes, we selected the substitution sequence contexts that were most representative of a given mutagenic process and examined their presence in the mutated genes. The APOBEC (APOBEC-like) mutational process was represented by TCW, (where W = A or T), and the 5-methylcytosine spontaneous deamination was represented by XCG (X = either of A, T, C, G). As expected, mutations representing the APOBEC signature were more frequent in the APOBEC cluster while the 5mC deamin-like type mutations were also higher in APOBEC and 5mC deamin cluster. Furthermore, the C>A mutations (dominant in SBS4) were enriched in the tobacco smoking cluster and the T>C mutations in the SBS16-like cluster. Nearly identical distribution was observed for all mutated genes as well as the subset of genes under selective pressure (Figure 5c), suggesting that the overall processes mutagenizing single bases identified for each cluster were also affecting the genes selected for during cancer development. Furthermore, we observed that the selected genes were relatively more enriched for small insertions/deletions (indels) compared to all mutated genes (Figure 5c). Collectively, these data indicate the key endogenous mutational processes acting primarily on cytosine residues that characterize the overall as well as driver gene mutation landscapes in NIRF oral cancers.

NIRF-enriched clusters exhibit keratinization and antimicrobial response gene expression programs

The causes of the SBS1 and SBS2/13 activity increase in the NIRF-enriched clusters are not known. To address this and to gain deeper functional and biological insights into the HNSCC patient stratification based on mutational signatures and the robustness of this approach, we performed DNA methylation analysis and differential gene expression and pathway analysis of the available TCGA methylome and transcriptomic data.

Previous studies have described distinct methylation profile in HNSCC subsites, including Brennan et al. who have reported five DNA methylation subtypes (Non-CIMP/CIMP Atypical, NSD1/ Stem-like Smoking and HPV+) using DNA methylation associated to mRNA expression (MethylMix). In addition to characterizing different epigenetic modifiers under selective pressure in the mutational signature-based clusters (Figure 5), we studied the distribution of the above methylation patterns in each cluster. Most patients in the two NIRF-enriched clusters harbored atypical, non-smoker DNA methylome patterns (Supplementary Figure S10).

Gene expression analysis of the NIRF-enriched clusters (see Methods) revealed higher expression of genes involved in the keratinization pathways (skin/epidermis development, keratinization, keratinocyte differentiation) (Figure 6a and b). In contrast, such activation of the keratinization pathway was absent in the tobacco smoking cluster (Supplemental Figure S11). Additionally, the NIRF-enriched clusters characteristically exhibited increased antimicrobial humoral response pathway activation (involving the genes S100A7, KLK7, CXCL11, CST9, WFDC12, and DEFB4A) (Figure 6a and b).

• Download figure
• Open in new tab

Figure 6:

Enrichment of the clusters for functional biological processes. (a) Pathway enrichment for genes differentially expressed in the 5mC deamination cluster (a) or APOBEC cluster (b) using the SBS16-like cluster as reference. Pathways with enrichment p-values < 0.05 are listed on the left alongside count indicators representing the number of genes enriched in each pathway. Network diagram (right) shows selected top pathways of interest and their respective differentially upregulated.

We next investigated whether the activation of APOBEC enzyme activity revealed by the APOBEC-directed mutagenesis is associated with a specific immune status of patients in the NIRF-enriched clusters. APOBEC-induced mutations in cancer have been shown to be generated by the APOBEC3A and APOBEC3B enzyme, with distinct isoform-specific substrate preferences ³¹. APOBEC3A preferentially targets the YTCA sequence (Y = pyrimidine bases) whereas APOBEC3B is more prone to mutate the RTCA sequences (R = purine bases). We analyzed the APOBEC mutation context and observed a higher APOBEC3A activity marked by the enriched YTCA context in the APOBEC cluster (Figure 7a). The APOBEC3A mRNA expression was elevated in the APOBEC cluster indicating that APOBEC3A remained active at the time of tumor sampling (Figure 7b). Interestingly, the APOBEC3A expression in the APOBEC cluster correlated with activated genes and pathways involved in the bacteria-specific immune response (Figure 7c). In sum, the mutational signature-based patient clustering coupled with gene regulation/expression analysis suggest a possible etiology of the NIRF OCSCCs associated with cell keratinization and immune responses against bacterial insults.

• Download figure
• Open in new tab

Figure 7:

The origin of APOBEC mutagenesis. (a) mutation counts in the APOBEC cluster for APOBEC3A (targeting the YTCA context, where Y stands for pyrimidine bases) and APOBEC3B (RTCA context, R standing for purine bases). (b) The APOBEC3A and APOBEC3B mRNA expression across clusters. (c) Top 20 pathways enriched for genes with levels correlated with APOBEC3A expression.

Samples and clinical data

United States HNSCC cancer incidence data were collected from Global Cancer Observatory (GCO) and Surveillance, Epidemiology, and End Results (SEER) websites [2022/12/02]. TCGA Clinical data (528 patients) were downloaded from the Genomic Data Commons [GDC release v33.0 downloaded on 2022/05/11]. Clinical, demographic and lifestyle variables (cancer subsite, age, sex, tobacco smoking, alcohol consumption) for all available whole exome sequenced (WES) samples were extracted. After exclusion of files corresponding to hypopharynx and oropharynx the dataset included OC (n=314) and LX (n=111) SCC samples. Clinical and demographic variables extracted from TCGA were as described in Table 1. Cases were classified into three age groups (young ≤ 40 years, middle-aged 41-69 years and old ≥70 years). Tobacco-smoking status was classified as lifelong non-smoker (less than 100 cigarettes in the lifetime) or ever smoker. Alcohol consumption status was classified as non-drinker (less than 1 drink/day) or ever-drinker. HPV infection status was obtained through the study of Campbell et al. ⁶³. Briefly, HPV status was assessed independently for all samples in two different centers by two different techniques (DNA sequencing and PathSeq algorithm and RNA-seq expression levels. Cases were classified as NIRF if they were non-smokers, non-drinkers, and HPV-negative.

Data pre-processing

Controlled-access single nucleotide variant files of 425 whole exome sequenced cancers from the TCGA-HNSCC cohort analyzed in this study were obtained from Genomic Data Commons [GDC release v33.0 downloaded on 2022/05/11], for each variant caller used available (MuSE, MuTect2, VarScan2 and Pindel). Single nucleotide polymorphisms and reads that did not pass the variant caller’s filters were removed. Somatic mutation spectra of samples were generated using SigprofilerMatrixGenerator (v1.2) with its R wrapper, SigProfilerMatrixGeneratorR (v1.1)¹³. Somatic mutations included single base substitutions (SBS), small insertions and deletions (indels or IDs).

Identification of rare mutational spectra

Signature bleeding may occur in the process of signature extraction when signatures identified in a small, highly mutated sample subset are ascribed to all samples. In order to detect mutation spectra outliers and mitigate potential signature bleeding, we first identified samples with high mutation burden by looking at mutation count distribution (Supplementary Figure S2a). Next SBS96 mutation spectra of samples with more than 600 mutations were compared using cosine similarity. Hierarchical clustering of those samples (using 1-cosine as distance) showed presence of nine distinct mutation spectra (Supplementary Figure S2b). Six of these spectra, with a mean spectrum exhibiting potential matches with COSMIC SBS signatures, were found in fewer than 5 samples and thus categorized as rare (Supplementary Figure S2c). One spectrum, identified in two lip cancers sample and one OC cancer sample, exhibited similarities to the UV light signature (SBS7). Additional two spectra (present in three samples) resembled various combinations of SBS signatures of DNA mismatch repair deficiency (SBS15, SBS20 and SBS21). One spectrum (in one sample) manifested with a strong C>G pattern enrichment similar to the COSMIC SBS39 signature. Finally, one spectrum (in one sample) showed no resemblance to a COSMIC counterpart, and another (in one sample) was only observed as a result of MuTect2 variant calling and is thus a possible artifact.

Extraction of de novo mutational signatures and sample clustering

After removing the atypical samples with high mutation count (see above), signature extraction was conducted on the 416 remaining samples using SigProfilerExtractorR (default parameters; v1.1.4) for each variant caller (MuSE, MuTect2, VarScan2) separately. While up to 10 de novo signatures were extracted in both the SBS96 and SBS384 catalogs, the optimal number of mutational signatures was assessed based on both the mean sample cosine and average stability, and set to 4 de novo signatures. Samples were finally clustered using K-means clustering algorithm into groups according to their signature relative contribution. Using the 384-channel signature dimension increased the stability across variant callers (Supplementary Figure S4a), by reducing the number of patients switching between clusters, from 21% (88/416) to 16% (67/416). Consequently, we consistently used the 384 channel-based signature analysis throughout the study.

A heatmap of SBS384 signature extraction was generated for each variant callers (Supplementary Figure S4b) and sample content in each cluster was compared among the variant caller-specific results, using Venn diagram (Supplementary Figure S4c). For each cluster, samples at the intersection of the 3 different variants were considered to be the most stable and were used for the final clustering. This final sample (n=347) assignment to each cluster was used in the remaining analysis steps.

COSMIC signature fitting

The all-COSMIC signature fitting by non-negative least squares (NNLS) was performed using MuTect2 variant caller-detected mutations, in samples with consistent clustering performance (see above). First, all COSMIC signatures were fitted to de novo extracted signatures using SigProfilerExtractorR (v1.1.4), using default parameters. Next, the MSA tool (v2.0) ²⁷ was used to fit the COSMIC signature set proposed by SigProfilerExtractor or to the well curated catalogue of known COSMIC signatures (version 3.2) to each sample individually. As no SBS384 version of COSMIC signature was available at the time of the analysis, the SBS288 transcription strand bias versions of COSMIC signatures was used in each case, with optimal threshold suggested by the MSA tool. Signature relative contribution was used to assess clustering stability. For final signature attribution, signature activities that included zero value inside the 95% confidence intervals were pruned to prevent overfitting.

Mutated genes under selective pressure

Genes under selective pressure were identified by computing non-synonymous to synonymous mutations ratio using dNdScv tool (v0.0.1.0) without covariate ⁶⁷. A dedicated mutation set was created for this analysis using Pindel in addition to the other variant callers used (MuSE, MuTect2, VarScan2). Only SBS mutations or indels identified by at least two variant callers were retained. A mean VAF value was calculated for each resulting mutation. TCW (with W standing for A or T) was used as template to identify APOBEC-like mutations and XCG (with X standing for A,T,C,G) for 5mC deamin-like mutation.

Differential Gene expression and methylation analysis

For gene expression analysis, open access RNA-Seq raw counts were download from GDC [GDC release v33.0 downloaded on 2022/05/11]. The analysis was focused on primary tumors exhibiting consistent clustering performance, resulting in 342 samples (see Extraction of de novo mutational signatures and sample clustering section).

Normalization and differential gene expression analysis of NIRF-enriched clusters (5mc deamin and APOBEC) were carried out against the other OC SBS16-like cluster using DESeq2 (version 1.38.3). As a control, the comparison of the tobacco smoking cluster to the SBS16-like cluster was performed.

Gene Ontology (GO) term enrichment analysis was conducted using the clusterProfiler package ⁶⁸ (version 4.6.2), on differentially expressed genes (DEG) with adjusted p-values less than 0.05 and a log fold change greater than 1.

Methylmix subtype were derived from conducted by Brennan et all. ⁶². The NSD1 mutation status was determined using Ensembl VEP annotation extracted from single nucleotide variant files. Only mutations with a high impact score were considered as inactivating mutations.

Role of funder

The funder had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Disclosures

The authors have no relevant information to disclose.