Significant efforts have been made over the years to better characterize and partition tumors into biologically and molecularly meaningful distinct subtypes. The rationale behind such tumor classification was to enable more precise, effective, and targeted therapeutic strategies. Such efforts have been rewarded, as in the case for breast cancer where death rates have dropped by 39% since 1989 (1), in part due to accurate subtype classification leading to efficacious treatments (2). This is particularly the case for the highly heterogenous group of triple negative breast cancers, where better understanding of the complex tumor microenvironment has facilitated tailored treatment regimens, including effective immunomodulation therapies (3, 4).

Similar subtype classification has been attempted for other tumor types with varied success. For example, head and neck cancers (HNSC) are clinically defined using the tumor-node-metastases (TNM) classification system and more recently the AJCC/UICC staging system, which rely heavily on pathological features such as primary tumor characteristics, lymph node spread, and distant metastasis, as well as relevant clinical history like smoking and alcohol usage (5, 6). Recent inclusion of using p16^INK4A status for HPV+ Oropharyngeal cancer in the eighth edition of TNM guidelines (7, 8) and PDL-1 expression paired with tumor mutational burden prior to treatment with pembrolizumab (9) represents examples where biomarkers have proven to be valuable for better diagnosis and therapy. Despite these advances, for many of the patients with advanced HNSC, histological and clinical staging do not correlate with clinical responses or prognosis (10). Therefore, continued efforts in subtype classification with the goal of more precise targeted therapeutics based on molecular signatures are needed.

Data mining of the TCGA datasets, primarily those based on transcriptomic signatures has resulted in subtyping of carcinomas from various anatomical sites. In the case of HNSC, tumors have been classified into four primary groups, Atypical (HPV+), Classical, Mesenchymal, and Basal subtypes (11, 12). However, given that HNSCs are inherently diverse and complex diseases of profound inter- and intra-tumoral heterogeneity, it is likely that additional subtypes exist. Support for such molecular diversity of HNSC comes from the discovery of a distinct subtype with impaired H3K36 function based on DNA methylation states (13) and from a meta-analysis approach that integrated multiple datasets (14). In general, classification of HNSC (and similarly other tumors) has by and large relied on a single facet of tumor biology, such as global gene expression profile, methylation patterns, somatic mutation states, or HPV infection status (15). Although such single-data driven approach has improved identification and our overall understanding of tumor subtypes, it is likely that application of an integrated bioinformatics-driven method that takes into account multiple data modalities might offer a powerful tool for new discoveries (16, 17).

Here, we establish a generalizable approach to integrate multimodal datasets for tumor subtype classification. Using HNSC as an example, we have performed a de novo integrated analysis of the bulk gene expression, methylation array, and microbiome datasets from TCGA Squamous Carcinoma patients. Specifically, we utilized multimodal spectral clustering followed by uniform manifold approximation and projection (UMAP), and potential of heat diffusion for affinity-based transition embedding (PHATE) dimension reduction (18–20). Our results demonstrate that data integration, followed by two-dimensional projection of the integrated similarity matrix perform robustly in partitioning samples, while preserving biological significance. Importantly, we show that during this process, the full range of molecular heterogeneity of HNSC tumors are captured, avoiding information loss. Furthermore, PHATE dimension reduction captures a dynamic gradient of expression in HNSC, ordering samples along a cellular trajectory toward more invasive Squamous Cell Carcinoma (SCC) as evident by tell-tale gene expression profiles. These results have broad implications for the challenge of accurate cancer subtype prediction by providing a robust strategy of multimodal data integration, leading to more accurate subtyping, with the potential for guiding future therapeutic intervention.

Large scale and integrated analysis of tumors, often based on multiple data modalities has revealed immense heterogeneity for most cancers. What has become quite clear from such studies is that cancers originating from the same human tissue or organ can dramatically differ in their etiology, pathology and the underlying genomic and epigenomic alterations. Precise tumor classification and identification of cancer subtypes based on underlying molecular signatures is thus an important step towards better understanding of the disease and importantly, for targeted therapy. This has led to recent studies such as the PanCancer Atlas integrative analysis, which identified 28 distinct molecular subtypes arising from the 33 different tumor types (21). Indeed, the unsupervised consensus clustering of tumor profiles performed across multiple TCGA genome-wide platforms including DNA methylation and RNA sequencing in the aforementioned study have reaffirmed the striking diversity of tumors and offered new insights into their molecular taxonomy. However, despite these advances, it is likely that there exist fine-grained tumor subtypes that are yet to be discovered. Such discovery would benefit from integrated application of new analytical and visualization tools to better parse the diverse sets of multi-dimensional data generated by high-throughput technologies.

Towards this end, several powerful strategies have been developed to interrogate multi-omic datasets (59). Examples of such strategies include low rank approximation based multi-omics data clustering (LRAcluster), which relied on linearly concatenate multi-omic profiles prior to probabilistic modeling and clustering of a latent subspace (60). This dimension reduction strategy allowed for integrative clustering of large-scale cancer multi-omics data. Other examples, including perturbation clustering for data integration and disease subtyping (PINS) and cluster-of-cluster assignments (COCA), which cluster each data view separately before they are integrated through consensus clustering strategies (61, 62). Similarly, statistical approaches using Bayesian statistics like canonical correlation analysis (CCA), iClusterBayes, BCC, are also capable of multi-modal analysis, with the added benefit of providing the probability of a sample belonging to a given cluster (63–66).

As an alternative to such integration algorithms, similarity-based strategies rely on first computing sample-to-sample measures of similarity or differences, prior to integration. This initial computation lends itself well to omic-based datasets of variable distributions (59). Recently, MoGCN has parlayed the advantage of autoencoders paired with neural networks for cancer subtype classification and analysis of breast invasive carcinoma (BRCA) samples (67). Indeed, similarity-based integration paired with spectral clustering techniques have provided an efficient method to handle complex genomic datasets of varying distributions (25, 26). Building upon these studies, here we present a computationally efficient analysis framework using similarity-based integration that allows for inter-cluster inspection and visualization at a much granular resolution.

We showcase the benefit of an integrated approach that relies on multimodality spectral clustering paired with dimension reduction techniques on large and complementary cancer genomic datasets. Indeed, we demonstrate the flexibility and performance of these algorithms, which in a seamless and robust fashion can handle datasets of varying size such as the pan-Cancer, pan-SCC and the HNSC patient cohorts from TCGA. Additionally, our strategy overcomes the limitation of traditional clustering approaches, which operate on single-view representations of a given sample set, and are designed to assign community participation based on a distance heuristic (68). This is particularly relevant for high dimensional data, like those analyzed in this study. Our results also suggest that the traditional classification approach of organizing tumor samples into discrete bins prior to calculating an average tumor profile based on a single mode of data offers only a static and limited view. Cancer progression is intimately associated with dynamic molecular changes in the molecular gene expression, DNA methylation states, and other modalities that are often transient in nature. We show that powerful dimension reduction techniques such as PHATE offer a better understanding of these transitions by ordering individual samples along these trajectories thus generating a non-discrete representation of the heterogeneity that is deeply inherent to most tumors.

We sought to probe the efficacy of dimension reduction and multi-view clustering strategies on HNSC tumors given that these represent distinct anatomical locations, diverse molecular mechanisms of carcinogenesis and highly heterogeneous tumor microenvironment. Unsurprisingly, patients with HNSC tumors have varying clinical responses among the 3-5 molecularly defined subtypes that have been primarily clustered anchored on single-data-view representations (69, 70). For our approach, we utilized integrated methylome, transcriptome, and microbiome signatures on TCGA HNSC tumors using recent advances in multi-graph integration using three data modalities. The nine clusters that were identified in our analysis included known HNSC subtypes such as the HPV+ and NSD1/H3K36 impaired groups, aided by inclusion of the tumor-specific microbiome, and DNA-methylation status, respectively. However, unlike previous reports, we found the basal subtype to be quite complex with four well-partitioned clusters that reflected a gradient of gene expression typified on one end by keratin enriched (B1) and mesenchymal markers (B4) on the other. These clusters also exhibited gradients in DNA methylation states and associated microbiome (Lactobacillus, Wolbachia, and Mannheimia among others) levels, further highlighting the power of the multimodal analysis. The incorporation of the microbiome state in our analysis is particularly interesting and worth further investigation because it is becoming clear that oral and gut microbiome are associated with HNSC development, progression, treatment and its potential side effects (71).

Several additional observations from our analysis of the HNSC tumors are note-worthy. First, by preserving global sample-to-sample relationships while ordering tumors by local affinities, PHATE highlighted the dynamic range of expression along the Epithelial-Mesenchymal transition (EMT) state in HNSC tumors. Furthermore. our study reaffirmed the evolving notion that tumor ecosystems often exhibit a continuum of meta-stable, intermediary pEMT (or hybrid) states between the epithelial and mesenchymal poles (72, 73). Second, clear segmentation of tumor samples characterized by heterogeneous cell populations and low tumor purity, such as stromal and immune cells, demonstrate the robustness and sensitivity of dimension reduction algorithms. We posit that unlike in the case of principal component analysis where the first two to three dimensions are driven by global variance, PHATE allowed us to model a broader and accurate representation of the multi-cellular tumor ecosystem. Indeed, this information is often lost using commonly used hierarchical clustering strategies which rely on pre-selected gene lists or only subsets identified a priori by principal component analysis.

Current treatment options for HNSC primarily consist of surgery, radiation and chemotherapy, administered in single or multi-modality regimens. However, these treatments still leave room for substantial improvement in efficacy and importantly does not consider molecularly defined subsets. Furthermore, approval of immunotherapy drugs such pembrolizumab and nivolumab for adjuvant treatment of recurrent and metastatic, or advanced unresectable HNSC, underscore the urgent need for research into tumor heterogeneity to better identify patients that are likely to be responders (74, 75). Given the complex link between immunomodulation and tumor microenvironment, a more robust and accurate classification strategy will be beneficial (76). Although we used HNSC to test our multimodal strategy, our studies provide a generalizable framework that can be applied to any cancer or other diseases. In addition, our approach can easily incorporate additional genomics data beyond methylation, transcriptome, and the microbiome. Finally, it is also important to consider novel data imputation strategies to maintain robust sample numbers in the event that not all data modalities are collected and available for each patient. Evolving techniques such as generative adversarial neural networks, among other machine learning strategies, can provide accurate predictions of both transcriptomic and methylation states that could be incorporated for multimodal clustering (77–79). Development of rapid and scalable multiomics data integration and mining strategies will continue to enable better modeling of the inherent heterogeneity of tumors and offer molecular insights into their complex and granular landscape.

Publicly available datasets were analyzed in this study. This data can be found here: https://www.cbioportal.org/comparison/overlap?comparisonId=60804602e4b0242bd5d4984c.

JB conceptualized the project, acquired and analyzed data, prepared the figures, and wrote the manuscript. NN provided funding and participated in the manuscript review and editing. MB and SS supervised the project, wrote and edited the manuscript, and provided funding. All authors contributed to the article and approved the submitted version.

This project was supported by grants from the National Institutes of Health (R01GM132199 to MJB, and R01AR073226 to SS) and Community Foundation for Greater Buffalo to MB, SS, and NN.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.