Highlights
- •Transcriptomic data can be transformed to 2D-images representing molecular pathways.
- •DeepClassPathway uses transcriptome 2D-images for explainable prediction.
- •DeepClassPathway allows explainable prediction of HPV-status in HNSCC.
Abstract
Objective
HPV-associated head and neck cancer is correlated with favorable prognosis; however, its underlying biology is not fully understood. We propose an explainable convolutional neural network (CNN) classifier, DeepClassPathway, that predicts HPV-status and allows patient-specific identification of molecular pathways driving classifier decisions.
Methods
The CNN was trained to classify HPV-status on transcriptome data from 264 (13% HPV-positive) and tested on 85 (25% HPV-positive) head and neck squamous carcinoma patients after transformation into 2D-treemaps representing molecular pathways. Grad-CAM saliency was used to quantify pathways contribution to individual CNN decisions. Model stability was assessed by shuffling pathways within 2D-images.
Results
The classification performance of the CNN-ensembles achieved ROC-AUC/PR-AUC of 0.96/0.90 for all treemap variants. Quantification of the averaged pathway saliency heatmaps consistently identified KRAS, spermatogenesis, bile acid metabolism, and inflammation signaling pathways as the four most informative for classifying HPV-positive patients and MYC targets, epithelial–mesenchymal transition, and protein secretion pathways for HPV-negative patients.
Conclusion
We have developed and applied an explainable CNN classification approach to transcriptome data from an oncology cohort with typical sample size that allows classification while accounting for the importance of molecular pathways in individual-level decisions.
Keywords
1. Introduction
One of the main applications of cancer omics is the prediction of prognostically relevant disease sub-entities. Many transcriptome-based classifiers predicting clinical outcome or defining transcriptional subtypes have been published in recent years, with most studies using a regression-based machine learning methodology [
[1]
]. Regression models come with the advantage of low complexity, which can be beneficial, but limits biological explainability. Deep learning with neural networks can make use of the full complexity of a dataset and allows reconstruction of the information used in decision making with explainable artificial intelligence methods.Strategies for using omics- and gene expression data as input to neural networks range from using them directly as gene expression vectors, representation vectors resulting from encoding/decoding or graph convolutional neural networks after mapping gene expressions to protein–protein interaction (PPi) networks [
[2]
]. However, these strategies do not sufficiently allow post hoc reconstruction of the biological information important for neural network decision making. One promising approach for overcoming this is to organize the gene expression into 2D-images with a structured representation of molecular pathways. Lopez-Garcia et al. (2020) first applied this concept by organizing gene expression into pathway 2D-treemaps as input for convolutional neural network (CNNs) predicting clinical outcome in a pan-cancer study [[3]
]. The resulting model showed good performance, but explainability was not investigated.High-risk human papilloma virus (HPV)-associated and HPV-negative head and neck squamous cell carcinoma (HNSCC), primarily associated with smoking and alcohol consumption, are regarded as distinct cancer entities which are clearly characterized by different clinical and molecular properties [
4
, 5
, - O'Sullivan B.
- Huang S.H.
- Su J.
- Garden A.S.
- Sturgis E.M.
- Dahlstrom K.
- et al.
Development and validation of a staging system for HPV-related oropharyngeal cancer by the International Collaboration on Oropharyngeal cancer Network for Staging (ICON-S): a multicentre cohort study.
Lancet Oncol. 2016; 17: 440-451
6
, 7
]. HPV-associated tumors treated with (chemo)radiation therapy have a better clinical prognosis and exhibit different transcriptional profiles compared to HPV-negative tumors. As comprehensive knowledge on the underlying biology of HPV- associated and -negative HNSCC exists and HPV-status typing is well established and robust we considered transcriptome data from HPV-characterized HNSCC as ideal real-world data for the development of an explainable deep learning approach for the prediction of disease.We developed and present DeepClassPathway, a workflow that transforms transcriptome profiles into 2D-images representing the gene expressions of 50 molecular pathways as an input for training and application of a classification CNN-ensemble (see Fig. 1 for an overview). Our concept was applied to a whole RNA-seq dataset from an in-house HNSCC cohort and a HNSCC cohort from The Cancer Genome Atlas (TCGA) [
[1]
,[4]
,[8]
]. The contributions of each molecular pathway to decision making were reconstructed by querying the magnitude of grad-CAM saliency heatmaps at the locations of pathways within the gene expression 2D-image [- Maihoefer C.
- Schuttrumpf L.
- Macht C.
- Pflugradt U.
- Hess J.
- Schneider L.
- et al.
Postoperative (chemo) radiation in patients with squamous cell cancers of the head and neck - clinical results from the cohort of the clinical cooperation group "Personalized Radiotherapy in Head and Neck Cancer.
Radiat Oncol. 2018; 13: 123
[9]
].
Fig. 1(color): Workflow of CNN-based prediction of disease subentities using treemap-transformed gene expression data. The upper panel describes the prediction of the HPV-status by a CNN-model using gene expression data transformed into pathway treemaps, followed by performance assessment using ROC-AUC (bottom right) and explainability by grad-CAM saliency mapping to the pathways (bottom left).
We describe a reproducible workflow for generating explainable CNN classification models, which can be generally applied to any transcriptome dataset with binary ground truth labels and demonstrate its predictive capabilities and the plausibility of the reconstructed biological information used in the model decision process.
2. Methods
2.1 HNSCC patient cohorts
The TCGA (n = 277) and the LMU-KKG (n = 72, Ludwig-Maximilians-University of Munich, Clinical Cooperation Group ‘Personalized Radiotherapy in Head and Neck Cancer’) HNSCC cohorts were analyzed [
[1]
,[4]
,[8]
]. Clinical patient data from cBioPortal and the TCGA GDC data portal [- Maihoefer C.
- Schuttrumpf L.
- Macht C.
- Pflugradt U.
- Hess J.
- Schneider L.
- et al.
Postoperative (chemo) radiation in patients with squamous cell cancers of the head and neck - clinical results from the cohort of the clinical cooperation group "Personalized Radiotherapy in Head and Neck Cancer.
Radiat Oncol. 2018; 13: 123
[10]
,[11]
] were used. The LMU-KKG cohort included patients with HNSCC of the hypopharynx, oropharynx, or oral cavity who had undergone surgical resection followed by adjuvant (chemo)radiotherapy treatment at the Department of Radiation Oncology, LMU University Hospital, Germany, between 2008 and 2013 [[8]
]. This study on clinical and biological data was approved by the local ethics committee in Munich (EA 448-13 and 17-116) and was carried out in accordance with the Declaration of Helsinki. HPV status for TCGA specimens was from the clinical patient data [- Maihoefer C.
- Schuttrumpf L.
- Macht C.
- Pflugradt U.
- Hess J.
- Schneider L.
- et al.
Postoperative (chemo) radiation in patients with squamous cell cancers of the head and neck - clinical results from the cohort of the clinical cooperation group "Personalized Radiotherapy in Head and Neck Cancer.
Radiat Oncol. 2018; 13: 123
[4]
]. HPV status for LMU-KKG specimens was determined as previously reported [[12]
].The clinical and histopathological data for TCGA and LMU-KKG patients included in the training and testing sets are presented in SI Table 1.
2.2 Transcriptome data preprocessing
A detailed description of transcriptome data preprocessing can be found in the supplementary extended methods section. In brief, the LMU RNA-seq data raw sequencing reads were adapter-trimmed, aligned and counted per gene. Raw counts of the TCGA cohorts were downloaded from the GDC data portal. All read counts were imported into DESeq2 before calculation of variance stabilization transformed (vst) gene expressions. The resulting vst-value matrix was used as input for generating the 2D-treemap images.
2.3 Transformation of gene expressions into 2D-treemap images
In a first step the gene expression vector (vst) of each patient was assigned to the genes defining each of the 50 MsigDB pathways and the resulting data frame was converted to 500 × 500 pixels treemaps. Each gene within the rectangles gets a greyscale value according to its gene expression value (high expression: dark, low expression: light). The MsigDB Hallmark genes have been selected because they contain only ground-truth evaluated genes [
[13]
]. As all 50 Hallmark pathways are contained, each 2D-image represents the molecular biological makeup of each tumor. To evaluate the robustness of our approach we generated three different variants of treemaps, leading to three different permutations of the pathway locations in the treemaps and optimized the CNN models separately for each one of the treemap variants.2.4 Deep learning model
We implemented a 2D-CNN [
[14]
,[15]
] with three convolutional blocks and 3 × 3 kernels for binary classification of HPV status. A detailed description of the 2D-CNN architecture can be found in the supplementary extended methods section.2.5 Model optimization and details
We trained the model using cross-validation on the training dataset, which was built by randomly selecting about 75% of the patients from the combined cohorts, with the remainder used as testing set. The training patients were subsequently split to perform 3-fold cross validation. This procedure was repeated 10 times with random initialization, thus 30 different CNN models were trained and validated. The models were then applied to the testing data and by averaging the predicted probabilities over the models, we obtained for every patient a probability for HPV-positive and HPV-negative. All models were trained using the categorical cross entropy loss function and the Adam optimizer. Weight decay was used additionally to dropout as regularization technique [
[16]
]. The hyper-parameters found to perform best are shown in Table 1.- Kingma D.P.
- Ba J.
Adam: a method for stochastic optimization.
arXiv:14126980. 2017; https://doi.org/10.48550/arXiv.1412.6980
Table 1Best set of CNN hyper-parameters found using the first treemap variant. The same hyper-parameters were then used for the other two variants.
| Hyper-parameter | Best performing value |
|---|---|
| Learning rate | |
| Weight decay | |
| Dropout rate | 0.30 |
| Batch size | 32 |
2.6 Pathway saliency mapping
To determine which pathways of the treemaps are related to HPV status, we used the gradient-weighted class activation mapping (grad-CAM) method [[
[9]
],[18]
]. It uses the gradients of a predicted class probability with respect to the feature set of the last convolutional layer to produce a coarse saliency map. The saliency maps are subsequently up-sampled to the shape of the input treemap, as can be seen in the bottom left part of Fig. 1. From the up-sampled saliency map, the saliencies of different tumor pathways were calculated by taking the average of the saliency map within all given pathway regions. Saliency maps were calculated separately for patients with clinically proven HPV-positive status with respect to the HPV-positive class and for patients with clinically proven HPV-negative status with respect to the HPV-negative class prediction. We computed pathway saliency only for models that predicted the correct HPV-status with scores higher than 0.70. Finally, the average saliency per pathway obtained for correctly classified patients was averaged over all models. Therefore, we obtained a patient and model averaged saliency per pathway and treemap variant (see Fig. 2).- Adebayo J.
- Gilmer J.
- Muelly M.
- Goodfellow I.
- Hardt M.
- Kim B.
Sanity checks for saliency maps.
arXiv:1810.03292. 2018; https://doi.org/10.48550/arXiv.1810.03292
Fig. 2(color): Mean grad-CAM saliencies obtained from the classification of HPV-positive (top) and HPV-negative (bottom) testing set patients. Shown are boxplots of grad-CAM saliencies averaged across the CNN model ensemble 1, sorted by median order from left to right.
2.7 Explainability
To visualize the relatedness of MsigDB Hallmark paths, a network graph was drawn, where pathways were represented as nodes and edges were weighted according to the number of overlapping genes between two paths. The Kamada Kawai layout was chosen so that the distance between nodes reflected their relatedness in terms of overlapping genes.
2.8 Predictive power
The discriminative performance in predicting HPV-status was evaluated with the receiver operating characteristics area under the curve (ROC-AUC) [
[19]
], and with the precision-recall area under the curve (PR-AUC) for imbalanced datasets [[20]
].3. Results
3.1 RNA-seq data
After random splitting, the training set (n = 264) contained n = 232 cases from the TCGA and n = 32 cases from the in-house LMU-KKG cohort and the test set (n = 85) consisted of n = 45 patients from the TCGA and n = 40 from the in-house cohort (SI Table 1). In total the dataset contained n = 22,866 genes.
3.2 Explainability
Heatmaps of the grad-CAM saliencies for each patient and model contained in model ensemble 1 can be found in SI File 1 for the prediction of HPV-associated and in SI File 2 for HPV-negative predictions. The top five mean grad-CAM salient MsigDB Hallmark pathways for the prediction of HPV-positive patients were spermatogenesis, KRAS-signaling (down), bile acid metabolism, inflammatory response and xenobiotic metabolism. In HPV-negative patients the most salient pathways were MYC targets (v1), apical surface, epithelial mesenchymal transition, MTORC1 signaling and reactive oxygen species (Fig. 2). The top salient pathways for both HPV-positive and negative for model ensemble 1, were also found in the top 10 of model ensemble 2 and 3 (SI File 3). The mean pathway saliencies in HPV-negative patients showed higher homogeneity expressed by smaller variation across pathways (SD: 0.017, range: 0.19–0.26) compared to HPV-positive patients (SD: 0.056, range: 0.08–0.34). The distributions of variances of the mean pathway saliencies in both groups were significantly different (Levene's homogeneity test p-value < 10−7).
The relatedness of the 50 Hallmark pathways in terms of genes shared between the underlying gene sets is visualized in Fig. 3, with mean grad CAM saliency values assigned to pathway nodes. The five major signaling pathways in HPV-positive predictions, KRAS-signaling (down), spermatogenesis, and inflammatory response pathways form a direct connection subnetwork, as do xenobiotic metabolism and bile acid metabolism. In HPV-negative predictions from the top five salient pathways MYC targets (v1), apical surface, epithelial mesenchymal transition, MTORC1 form a directly connected subnetwork with no direct connection to the reactive oxygen species pathway. In order to identify grad-CAM saliency patterns the mean grad-CAM saliencies per patient and pathway were subjected to unsupervised hierarchical clustering, and four main grad-CAM patterns were found (SI Fig. 2).
Fig. 3(color): Hallmark pathway gene overlap network with mapped mean grad-CAM saliencies of HPV-positive (upper panel) and HPV-negative (lower panel) classified testing set patients. Network nodes are Hallmark pathways colored by the grad-CAM saliencies averaged across CNN model ensemble 1 and patients. Edges indicate shared genes between two paths, and the edge length is weighted by the number of shared genes. The closer the network nodes are to each other, the more genes they have in common. The top five salient pathways are indicated by pink circles. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
3.3 Predictive power
For the three treemap variants and the corresponding ensembles of models, ROC-AUC was 95.5%, 95.9% and 95.8% and PR-AUC was 90.1%, 90.4% and 90.4% (Table 2 and SI Fig. 1).
Table 2Prediction performance of CNN models for different input treemap variants.
| CNN-Model | ROC-AUC test | PR-AUC test |
|---|---|---|
| Model ensemble 1 | 0.955 | 0.900 |
| Model ensemble 2 | 0.959 | 0.904 |
| Model ensemble 3 | 0.958 | 0.904 |
4. Discussion
In addition to assessing predictive power, which was consistently good across model ensembles, we thoroughly explored the capabilities of our approach with respect to explainability expressed by grad-CAM. Overall, we observed higher homogeneity of pathway mean grad-CAM saliencies in HPV-negative- compared to HPV-positive patients. This goes along with a clear ranking of mean pathway grad-CAM saliencies observed in HPV-positive- but not in HPV-negative-predicted patients. In contrast to HPV-negative HNSCC, tumorigenesis of HPV-related tumors is characterized by deregulated expression of the HPV oncogenes E6 and E7, which inactivate RB and p53, respectively, and lead to immortalization and cancer progression through the accumulation of typical alterations such as PIK3CA mutations or amplifications [
[6]
,[7]
]. The process of tumorigenesis in HPV-negative tumors is not triggered by HPV viruses and is therefore rather stochastic, which on average should lead to a more regular distribution of activated pathways as observed in our data. The top five pathway grad-CAM saliencies in HPV-associated predicted patients were spermatogenesis, KRAS signaling, bile acid metabolism, inflammatory response and xenobiotic metabolism. Martinez et al. (2007) reported specific upregulation of spermatogenesis genes in HPV-positive HNSCC [[21]
]. Further, it is known that HPV-driven HNSCC regulate the inflammatory response in the frame of immune escape [[22]
]. Only little knowledge is available on disturbed KRAS-signaling in HNSCC—only rare occurrence of KRAS-mutation is reported [[23]
]. Bile acid metabolism and xenobiotic metabolism are not plausible in the context of HNSCC. This can be explained by the unspecific detection of their saliencies due to genes present not only in these but also in other signaling. Hence, improving the specificity of ground truth-based pathway gene sets would be an important next step towards improving their usability. The pathways MYC targets (v1), apical surface, epithelial mesenchymal transition, MTORC1 signaling and reactive oxygen species with highest grad-CAM saliencies in the HPV-negative predicted patients are consistent with published findings on HNSCC and squamous cell carcinogenesis [23
, 24
, 25
, - Kim S.W.
- Cheon K.
- Kim C.H.
- Yoon J.H.
- Hawke D.H.
- Kobayashi R.
- et al.
Proteomics-based identification of proteins secreted in apical surface fluid of squamous metaplastic human tracheobronchial epithelial cells cultured by three-dimensional organotypic air-liquid interface method.
Cancer Res. 2007; 67: 6565-6573
26
] and therefore highlight the plausibility of the approach. The pathways which showed highest mean grad-CAM activations showed overall consistent patterns in HPV-associated and HPV-negative patients as identified by hierarchical clustering. Two out of three HPV-associated patients misclassified as HPV-negative by our models showed an overall low grad-CAM activation which is probably due to faulty transcriptomic profiles. The third HPV-associated patient misclassified as HPV-negative and the HPV-associated cases classified as HPV-negative showed grad-CAM patterns which are typical for their predicted class, which could be due to clinical misclassification.5. Conclusion
In this study, we developed a deep learning model for the prediction of HPV-status from transcriptome data in head and neck cancer patients from two independent cohorts which achieved very good and stable classification performance. The advantage of the proposed CNN model was identified in its biological explainability, as it provides the molecular pathway information that determined the classification at an individual patient level. In conclusion, we showed that CNNs can be successfully applied to transcriptomic datasets with sample sizes usual in translational clinical science to predict a meaningful endpoint and obtain information about the biological pathways driving the model's prediction. This provides the starting point for the further development of explainable transcriptome-based CNN prediction models.
Data availability
The TCGA HNSC data used in this paper are available from the GDC data legacy portal. RNA-seq data generated from the in-house HNSCC cohort are available from Gene Expression Omnibus under the accession number GSE205308.
Funding
This study was supported by the German Federal Ministry of Education and Research ( BMBF ) project ZiSStrans ( 02NUK047A to JH, KU and HZ and 02NUK047C to CB and KL).
Author contributions
Conceptualization, G.L., J.H., and K.U.; methodology, C.K., E.L., G.L., J.H., M.R. and K.U.; software, E.L. and K.U.; validation, C.K., J.H, E.L., G.L., M.R. and K.U.; formal analysis, E.L., J.H., and K.U.; investigation, E.L., J.H., K.U., S.M. and U.P.; resources, A.W., C.B., F.K., H.Z., M.C. and P.M.; data curation, E.L., J.H., S.M., K.U. and U.P.; writing—original draft preparation, E.L., G.L., J.H. and K.U.; writing—review and editing, C.B., C.K., F.K., G.L., K.L., M.R. and S.M.; visualization, E.L., J.H. and K.U.; supervision, C.B., C.K., G.L., M.R. and K.U.; project administration, G.L. and K.U.; funding acquisition, C.B., H.Z. All authors have read and agreed to the published version of the manuscript.
Conflict of interest statement
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
The authors want to thank S. Heuer, L. Holler and C. Innerlohinger for their excellent technical assistance; all co-workers of the Clinical Cooperation Group “Personalized Radiotherapy in Head and Neck Cancer” for scientific support.
Appendix A. Supplementary data
The following are the Supplementary data to this article:
SI Fig. 1: Model averaged ROC-AUC and PR-AUC HPV-status prediction for the CNN models obtained on the testing dataset. Model numbers indicate different input treemap variants.
SI Fig. 2: Grad-CAM saliencies averaged across CNN model ensemble. Pathways (rows) are sorted according to hierarchical clustering with four main clusters colored in red, blue, green and yellow and columns (patients) sorted after ground-truth HPV status (left: HPV-negative, right: HPV-positive). Ground-truth and CNN-predicted HPV-status and misclassified cases are indicated by color bars.
- Multimedia component 1
- Multimedia component 2
- Multimedia component 1
- Multimedia component 4
References
- A five-MicroRNA signature predicts survival and disease control of patients with head and neck cancer negative for HPV Infection.Clin Cancer Res. 2019; 25: 1505-1516
- Deep learning in cancer diagnosis, prognosis and treatment selection.Genome Med. 2021; 13: 152
- Transfer learning with convolutional neural networks for cancer survival prediction using gene-expression data.PLoS One. 2020; 15e0230536
- Comprehensive genomic characterization of head and neck squamous cell carcinomas.Nature. 2015; 517: 576-582
- Development and validation of a staging system for HPV-related oropharyngeal cancer by the International Collaboration on Oropharyngeal cancer Network for Staging (ICON-S): a multicentre cohort study.Lancet Oncol. 2016; 17: 440-451
- Clinical update on head and neck cancer: molecular biology and ongoing challenges.Cell Death Dis. 2019; 10: 540
- The molecular landscape of head and neck cancer.Nat Rev Cancer. 2018; 18: 269-282
- Postoperative (chemo) radiation in patients with squamous cell cancers of the head and neck - clinical results from the cohort of the clinical cooperation group "Personalized Radiotherapy in Head and Neck Cancer.Radiat Oncol. 2018; 13: 123
- Grad-CAM: visual explanations from deep networks via gradient-based localization.Int. J. Comput. Vision. 2019; 128: 336-359https://doi.org/10.1007/s11263-019-01228-7
- The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data.Cancer Discov. 2012; 2: 401-404
- Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal.Sci Signal. 2013; 6: pl1
- A prognostic mRNA expression signature of four 16q24.3 genes in radio(chemo)therapy-treated head and neck squamous cell carcinoma (HNSCC).Mol Oncol. 2018; 12: 2085-2101
- The Molecular Signatures Database (MSigDB) hallmark gene set collection.Cell Syst. 2015; 1: 417-425
- Deep learning.Nature. 2015; 521: 436-444
- ImageNet classification with deep convolutional neural networks.Adv Neural Inf Process Syst. 2012; 25
- Adam: a method for stochastic optimization.arXiv:14126980. 2017; https://doi.org/10.48550/arXiv.1412.6980
- Sanity checks for saliency maps.arXiv:1810.03292. 2018; https://doi.org/10.48550/arXiv.1810.03292
- The meaning and use of the area under a receiver operating characteristic (ROC) curve.Radiology. 1982; 143: 29-36
- The precision-recall plot is more informative than the ROC plot when evaluating binary classifiers on imbalanced datasets.PLoS One. 2015; 10e0118432
- Identification of differentially expressed genes in HPV-positive and HPV-negative oropharyngeal squamous cell carcinomas.Eur J Cancer. 2007; 43: 415-432
- The Double-Edged Sword-how human Papillomaviruses Interact with Immunity in head and neck cancer.Front Immunol. 2019; 10: 653
- Molecular characterization of head and neck cancer: how close to personalized targeted therapy?.Mol Diagn Ther. 2012; 16: 209-222
- Epithelial-to-mesenchymal transition-derived heterogeneity in head and neck squamous cell carcinomas.Cancers (Basel). 2021; 13
- Proteomics-based identification of proteins secreted in apical surface fluid of squamous metaplastic human tracheobronchial epithelial cells cultured by three-dimensional organotypic air-liquid interface method.Cancer Res. 2007; 67: 6565-6573
- Distinct subpopulations of head and neck cancer cells with different levels of intracellular reactive oxygen species exhibit diverse stemness, proliferation, and chemosensitivity.Cancer Res. 2014; 74: 6291-6305
Article Info
Publication History
Published online: September 30, 2022
Accepted:
August 25,
2022
Received:
July 7,
2022
Identification
Copyright
© 2022 Elsevier Ltd. All rights reserved.
00533-0&id=gr1.jpg){kind=link}
00533-0&id=gr2.jpg){kind=link}
00533-0&id=gr3.jpg){kind=link}