{'fraction_non_alphanumeric': 0.01739664976179499, 'fraction_numerical': 0.010388812048563087, 'mean_word_length': 5.271395528141865, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A disease network‐based deep learning approach for characterizing melanoma

Multiple types of genomic variations are present in cutaneous melanoma and some of the genomic features may have an impact on the prognosis of the disease. The access to genomics data via public repositories such as The Cancer Genome Atlas (TCGA) allows for a better understanding of melanoma at the molecular level, therefore making characterization of substantial heterogeneity in melanoma patients possible. Here, we proposed an approach that integrates genomics data, a disease network, and a deep learning model to classify melanoma patients for prognosis, assess the impact of genomic features on the classification and provide interpretation to the impactful features. We integrated genomics data into a melanoma network and applied an autoencoder model to identify subgroups in TCGA melanoma patients. The model utilizes communities identified in the network to effectively reduce the dimensionality of genomics data into a patient score profile. Based on the score profile, we identified three patient subtypes that show different survival times. Furthermore, we quantified and ranked the impact of genomic features on the patient score profile using a machine‐learning technique. Follow‐up analysis of the top‐ranking features provided us with the biological interpretation of them at both pathway and molecular levels, such as their mutation and interactome profiles in melanoma and their involvement in pathways associated with signaling transduction, immune system and cell cycle. Taken together, we demonstrated the ability of the approach to identify disease subgroups using a deep learning model that captures the most relevant information of genomics data in the melanoma network.

identification of genomic features most relevant to patient outcome. Here, the authors integrate genomics data with a disease network and deep learning model for the prognostic classification of melanoma patients and assessment of impacts of genomic features on disease classification.
The data suggest that deep learning models based on genomics data and disease networks can contribute to personalized prognostic assessment for melanoma patients. The generic nature of the approach suggests that it is applicable to other cancer types.

| INTRODUCTION
Cutaneous melanoma is one of the most aggressive forms of skin cancer and is characterized by a lymphogenic and hematogenic metastatic spread which is directly proportional to the vertical depth of invasion of the tumor cells into the skin. 1 Early detection and surgical excision of the primary tumor are essential as the disease becomes life threatening once it metastasizes. 2 The introduction of small-molecule inhibitors targeting the kinases BRAF and MEK as well as immune checkpoint inhibitors have tremendously changed the treatment landscape in the last 10 years, leading to remarkable treatment results in patients with metastatic melanoma. 3 Nevertheless, many patients with advanced melanoma still do not benefit from these new drugs either through primary or acquired treatment resistance. 4 Despite enormous research, the understanding of melanoma progression and certain melanoma subtypes remains very limited. The Cancer Genome Atlas (TCGA) database provides information on copy number, mutation, transcriptomic, methylation and clinical data of melanoma patients. These data allow integrative analysis of the genomic profiles of melanoma and thereby improves our understanding of the progression and development of the disease. Prominent studies include the integration of genomics data to identify novel driver genes of melanoma 5 and informative individual markers and pathways that provide valuable insights into melanoma prognosis. 6 Besides, tumor heterogeneity impedes the efficacy of available therapies for melanoma, 7 imposing additional needs for developing tools that can stratify melanoma patients that show clinical and genomic distinction, therefore facilitating the development of precision medicine.
Deep learning models, which allow efficient extraction of complex data with multiple levels of abstraction, have proven a powerful tool in facilitating the classification of cancer and improving cancer prognosis by integrating genomics data and clinical factors. 8 For instance, Xie et al fed four types of data into a model combining a multimodal neural network with an autoencoder and identified four genomic clusters in 8646 patients with 29 cancer types. 9 Applying the identified genomic clusters to melanoma cohorts treated by anti-  showed that patients with distinct genomic features have different responses to the immunotherapy. Specifically, the patient subtype with high microsatellite instability and somatic copy number alterations but low tumor mutation burden showed the lowest clinical benefit of the immunotherapy and the shortest survival. Although there are variant architectures in deep learning models, most of them consist of multiple fully connected layers for analyzing genomics data, making them difficult to provide biological interpretation. Hence, researchers added annotated information such as pathways into neural networks to increase the interpretability of deep neural networks.
Hao et al developed a pathway-based sparse deep neural network that integrates genomics and clinical data to identify prognostic factors for glioblastoma and serous cystadenocarcinoma. 10 The addition of a pathway layer between the input layer (ie, the genomics data) and the hidden layer allowed biological interpretation of the identified genomic features at the pathway level. Similarly, Feng et al built a deep neural network by connecting the transcriptomic and copy number variation profiles of 1967 genes with 46 cancer pathways to predict cancer patients' survival time. 11 Computation of scores of the pathways quantified their impact on the model outputs and thus explaining the association between the pathways and the patients' survival time. The interpretable deep learning models not only leverage the full potential of genomics data but also provide biological explanations that are crucial for understanding tumor heterogeneity, therefore allowing for identifying therapeutic targets for specific patient groups.
Additionally, carcinogenesis can be viewed as a rewired network due to endogenous perturbations resulting from genomic alterations.
Thus, researchers integrate human genomics and network analysis to refine our biological understanding of cancer and identify novel prognostic markers. One advantage of network analysis is that interactome information can be stratified into communities based on network topology, thus permitting a more precise understanding of molecular mechanisms. The success of such approaches has been demonstrated in predicting patient-specific mutation signatures for melanoma, 12 identifying synergistic drug combinations for diverse melanoma genomic subtypes, 13 and detecting cooperative microRNAs (miRNAs) for treating metastatic melanoma resistant to anti-PD-1 therapy. 14 Thus, it is appealing to develop an approach that integrates genomics data and a disease network with a deep learning model to characterize melanoma.
In our study, we integrated genomics data into a melanoma network and applied a deep learning model to identify subgroups in the TCGA melanoma patients (Figure 1). The model utilizes communities identified in the melanoma network to effectively reduce the dimensionality of genomics data into a patient score profile. Using the score profile, we identified three subtypes of melanoma patients. Post hoc analysis of the patients' survival time and clinical features provided supporting evidence and interpretation of the identified clusters. Furthermore, we applied a machine-learning technique to quantify and rank the impact of genomic features on the score profile. We integrated gene set enrichment analysis with the melanoma network to perform functional analysis that provides the biological interpretation of the top-ranking genomic features at the pathway and molecular levels. By doing so, we identified a group of genes (eg, MMP2, PPP3CA and TUBB2B) that are impactful on the patient score profile, characterized their involvement in signaling pathways and depicted their molecular features such as mutation and interactome profiles in melanoma. Taken together, the approach allows for the effective integration of genomics data into a disease network to capture data variation for robust identification of disease subgroups with distinct clinical and molecular features.

| The melanoma network and network community detection
The melanoma map was constructed in our previous study and downloaded from www.vcells.net/melanoma. 15 Briefly, the core map was created and curated through literature reading and annotation (see Supporting Information Materials for details). The principle of network community detection is to find a series of subnetworks with high cohesion and low coupling (see Supporting Information Materials for details). Based on this principle, we used a simulated annealing algorithm to partition the melanoma network. 16 As the algorithm can produce different network partitions, we ran numerous simulations on the melanoma network and ranked them based on network modularity scores. The top 10 partitioned networks have similar modularity scores and a comparable number of communities ( Figure S1), suggesting that the identified communities represent the optimal partition of the network. We selected the fully connected partitioned network with 21 communities and the maximum modularity score (0.485) for subsequent analysis. We ran the simulations and performed the analysis using the igraph package in R.

| TCGA data processing
We downloaded the TCGA skin cutaneous melanoma dataset from the

| Network community scores
We adapted a sparse autoencoder model called PathME to transform genomics data of melanoma patients into a network community-based F I G U R E 1 Workflow of the study. We integrated the TCGA data into a melanoma network and performed community identification in the network. The genomic profiles of network community genes were inputs of an autoencoder model that transforms the genomics data into a community score. The scoring profile of patients was used for patient classification followed by subsequent analysis including patients' survival analysis, ranking of genomic features and interpretation of the features [Color figure can be viewed at wileyonlinelibrary.com] score profile. The model has an unsupervised multimodal neural network architecture that can learn a low-dimensional embedding of genomic features from multiple sources. 17 The structure of the model is shown in Figure 3, and a detailed description of model configuration and training can be found in Supporting Information Materials. After training the model for the 21 identified network communities sequentially, we concatenated scores of network communities for each patient. This results in a community-based score matrix (21 rows and 463 columns) that represents a reduced representation of the input data and can be used for patient clustering.

| Patient stratification
We performed agglomerative hierarchical clustering using network community scores to identify clusters in patient samples. The algorithm first treated each patient sample as an individual cluster and calculated distances between clusters. The two samples with the shortest distance were merged into one cluster, and this process was repeated until all patient samples were merged into the same cluster.
For calculating distances between clusters that have more than one sample, we used the complete-linkage clustering method to calculate their distances. The optimal number of clusters was determined by 26 indices (such as silhouette and C-Index) that are well-established metrics to assess cluster reliability. This analysis was performed using Nbclust package in R.

| Shapley additive explanations
Shapley additive explanations (SHAP) is a game theoretic approach that provides interpretation to the output of machine learning models. 18 1): where

| Survival analysis and statistical tests on clinical features
We applied chi-squared analysis to test interdependences between the patient clusters and their clinical information such as gender, age, tumor types, tumor stages and whether receiving any treatment (ie, radiation and pharmaceutical therapy). The number of clustered patients in different categories of clinical features was used for this test. We applied Kruskal-Wallis tests to detect interdependences between the patient clusters and their genomic profiles such as somatic mutation frequency and copy number variation frequency.
The somatic mutation frequency of a patient is a sum of counts of any type of mutations identified in all genes in the melanoma network.
Copy number variation frequency of a patient is the number of loss and gain in copies of all genes in the melanoma network. The Kruskal-Wallis test was used for comparing all three clusters and any two of Note: The table shows the way of processing the data. The raw data were first mapped from identifiers to genes and then integrated into the melanoma network, resulting in four matrixes for further analysis.
the three clusters. P-value ≤.05 was considered significant. The statistical analysis was performed in R.
Patients' survival time after the diagnosis of melanoma was used for performing survival analysis. We performed log-rank tests to compare the survival curves among the clustered patients. In 453 of 463 patients, survival information was available in TCGA data and used to draw the Kaplan-Meier plot using the TCGAbiolinks package.
The mutation annotation format file of melanoma from TCGA was obtained using TCGAbiolinks. We draw somatic mutation profiles of genes of interest using maftools in R.

| Gene set enrichment analysis
To identify pathways for which genes in each network community are enriched, we performed gene set enrichment analysis using the Reactome database. 19 The pathways in the database were classified into 26 root categories, such as signaling transduction, immune system and cell cycle. 20 P-values of Fisher's exact test were corrected using the Benjamini-Hochberg method. Pathways with a false discovery rate (FDR) ≤ 0.05 were regarded as significant. This analysis was performed using ReactomePA package in R. In addition, we selected pathways containing at least 10 genes for drawing heat maps that show top-ranking features and the corresponding genomic information and drug-target interaction in different pathway categories. The heat maps were drawn using Complexheatmap in R.

| Data visualization
The melanoma network annotated with the identified communities and the networks for the top-ranking features were visualized using Cytoscape 3.82. The individual network communities annotated with genomic profiles were drawn using igraph and ggplot2 in R.

| Community identification in the melanoma network
We used a manually curated map that accounts for state-of-the-art knowledge on the molecular biology of melanoma and contains 16 signal transduction pathways, such as growth factor receptor signaling, melanocortin signaling, death-ligand signaling, TNFα signaling and P53 signaling. 15 We expanded the map by considering miRNA-gene interactions, interactions ( Figure 2A). The topology of the network is scale free, as its degree distribution approximates a power law ( Figure 2B).
We integrated the network with genomics data of 463 melanoma patients, including transcriptome, DNA methylation, somatic mutations and copy number variation (see Section 2). Such a network provided us with a comprehensive view of genomic profiles and molecular interactome underlying the pathogenesis and development of melanoma ( Figure 2C). By using a simulated annealing method, we identified 21 network communities with a modularity score of 0.485 (see Section 2; Figure 2). The size of the identified communities ranges from 6 to 579 nodes with the number of edge size ranging from 6 to 1914. Among the communities, 18 are scale free and the rest have a star topology ( Figure S3). The identified network communities annotated with melanoma patients' genomic profiles allowed us to cluster patients and analyze their clinical characteristics.

| Patient stratification based on communitybased scores
An autoencoder is a type of artificial neural network that can be used for dimensionality reduction and feature learning in an unsupervised manner. 21 Such a neural network is composed of two parts: an encoder that maps input data to a hidden layer and a decoder that maps the hidden layer to an output layer. An autoencoder can reconstruct data by minimizing the error between the input layer and the reconstructed layer. Specifically, if the size of the hidden layer is less than the size of the input layer, the hidden layer can capture the significant features of the input data.
To detect potential clusters in melanoma patients, we fed the identified network communities with gene expression and DNA methylation data into an autoencoder model (see Section 2; Figure 3). For the reason that the mutation and copy number variation data are sparse and not continuous, we did not consider them as inputs of the autoencoder model but used them for analyzing the genomic differences in clustered patients. The model reduced the dimension of the genomics data by transforming them into a network-community score profile, which is a matrix with 21 communities and 463 patients. Next, we applied hierarchical clustering on the community score profile to classify the patient samples and obtained three patient clusters (see Section 2). Such classification was systematically estimated and supported by 14 indices out of 26 ( Figure 4A). In addition, the identified clusters showed a clear pattern in their community scores that were not seen using genomics data ( Figures 4B and S4). Compared to the individual genomics data, the community score profile showed better classification of patient samples while applying a nonlinear dimensionality reduction method to embed the high-dimensional data in a 2-dimensional space ( Figure 4C). Taken together, the results suggested that the community score has the advantage in identifying subgroups in patients compared to using the genomics data alone. The identified patient clusters may have different clinical features that require further analysis.

| Association between the patient clusters and their clinical features and genomic profiles
We performed statistical tests to investigate whether the identified patient clusters have differences in clinical features and genomic profiles such as somatic mutation and copy number variation (see Section 2). The identified patient clusters showed significant differences in survival time-Cluster 3 has a much shorter median survival time and lower surviving probability than the other two clusters ( Figure 4D). Besides, the clusters showed a significant difference in tumor types-Cluster 3 has a lower percentage of patients with metastasis than the other clusters ( Table 2). No significant differences among the clusters were identified in other clinical features such as patients' gender, age, tumor stage and treatment differences (Table 2).
Furthermore, Clusters 1 and 2 showed significant differences in frequency of somatic mutation, and no significant difference was identified in copy number variation between the patient clusters ( Table 2; Figure S5). Taken together, the data showed that the identified patient clusters have differing survival times and different profiles in terms of tumor types and somatic mutation. Thus, it is interesting to explain such differences at the molecular level.

| Interpretation of genomic features learned for patient stratification
We combined the SHAP method and gene set enrichment analysis to interpret the genomic features learned by the autoencoder model for patient stratification. The former quantified the importance of genomic features (see Section 2), and the latter helps to explain the function of the network community genes at both pathway and molecular levels. For most network communities, the gene expression data show a much stronger contribution to the community score used for patient clustering than the methylation data, as the corresponding SHAP values are much higher ( Figure 5). The only exception is community 12, and its features in the methylation data have much higher SHAP values than the features in the gene expression data. Besides, the contribution of the molecular features on the community scores is similar for individual patient clusters ( Figure S6). (1 node) Patient classification The architecture of the autoencoder model. The model has one input layer, three hidden layers and one reconstructed layer. The input features are the z-score normalized gene expression and DNA methylation profiles of nodes in the identified network communities. The input layer is connected to the first hidden layer using tanh activation function. The second hidden layer has only one node that represents reduced dimensions of the input data for the network community. The decoder is symmetrical to the encoder where the node in the second hidden layer is split into two parts to form the third hidden layer. The third hidden layer is connected, via tanh activation function, to the reconstructed layer (ie, the output layer) that has the same dimension as the input layer. We trained the model sequentially using the genomics data of the 21 network communities and obtained a community-based score matrix (21 rows and 463 columns) that represents a reduced representation of the input data. We further used the score matrix for patient clustering. The data showed that the cost function converges to a minimum for each community ( Figure S2). The domain of hyperparameters and the best hyperparameter set for each network community are shown in Tables S2 and S3, Table S1 for details). For the genomics data, the table shows the corresponding P-value of the Kruskal-Wallis test (KWT) for the overall comparison between three clusters and individual comparisons among them.
0.00 0.05 0.10  Figure 7C). Besides, TUBB2B is involved in the Hedgehog signaling pathway whose activity plays a crucial role in melanoma tumorigenesis. 43 TUBB2B has a higher copy number variation frequency in the Cluster 3 patients than the other two clusters, and it shows mutations in only two melanoma patients that belong to Clusters 1 and 3, respectively ( Figure S8C). It has been shown that the mRNA and protein expression levels of TUBB2B are upregulated in metastatic melanoma cells, suggesting its important role in regulating melanoma. 44 Besides, the network showed that TUBB2B is targeted by miR-93a-3p and such post-transcriptional regulation is identified using a high throughput technique called CLASH that can recover RNA duplexes bound by Argonaute protein ( Figure S9B).  2 3  1 3  7  1 3  5  86 2  22  5  27  1 1  3  1  1 0  5  2  3  3  63  16  11  3 2  12  35 5  51  2  13  10  16  18  2 1  19 6  3  3 2  5  1  1  1  5  2 45 MMP16 is highly expressed in melanoma metastases but not in primary tumors 46 and has been identified as an effector that promotes melanoma cell invasion depending on the matrix composition. 46 Additionally, MMP16 overexpression has been linked to increased lymphatic invasion and a worse clinical prognosis of melanoma patients, 46 and it seems to cooperate with melanoma chondroitin sulfate proteoglycan (MCSP) to promote melanoma invasion. 47  recently. S100B ranked first in the methylation data of network community 13. This calcium-binding protein has been previously described as an independent prognostic biomarker in melanoma because its elevated baseline serum levels are associated with poor overall survival. 50 EIF1AX, which ranked sixth in gene expression data of network community 15, is frequently mutated in uveal melanoma, 51 a rare melanoma subtype that was not included in the analyzed dataset.
Mutations in EIF1AX lead to an aberrant translational regulation and have been associated with a low risk for developing metastasis in uveal melanoma. 52 ERBB3 ranked ninth in gene expression data of network community 16 and encodes Erb-B2 receptor tyrosine kinase 3, a frequently overexpressed member of the epidermal growth factor receptor (EGFR) family of receptor tyrosine kinases in melanoma. 53 ERBB3 is associated with adaptive drug resistance in melanoma as its activation mediates resistance to BRAF and MEK inhibitors by activating PI3K/AKT signaling. 54 Interestingly, its expression can be directly regulated by the transcription factor SOX10, which is also among the top-ranking features of network community 7. 55 ERBB3 inhibition using neutralizing antibodies enhances the efficacy of MEK inhibitors in BRAF and NRAS wild-type melanoma. 56 Additionally, ERBB3 can be targeted by the small molecule inhibitor sapitinib (AZD8931) that also inhibits ERBB1 and ERBB2. 57 The neural crest stem cell marker NGFR ranked first in gene expression data of network community 21.
The gene is upregulated in dedifferentiating melanoma cells that acquire a neural crest stem cell-like phenotype as a response to BRAF and MEK inhibition. 58 Taken together, combining the SHAP method, the melanoma network and the gene set enrichment analysis provides us with a powerful tool to interpret the impactful features that are learned by the autoencoder model to stratify melanoma patients. This not only elaborates the molecular function of the relevant genes but also unravels their association with signaling pathways that underlie the pathogenesis and development of melanoma. Furthermore, analyzing the genomic features of the patient clusters helps to explain differing survival times and potential differences in treatment response.

| DISCUSSION AND CONCLUSIONS
The rapid production and easy access of genomics data have created an unprecedented opportunity to study cancer biology. However, multidimensional data pose a huge challenge for data analysis. In our study, we developed an integrative approach that combines genomics data, a disease network, and a deep learning model to identify and analyze subtypes of melanoma patients. The genomics data of melanoma have high dimensionality, is noisy, and contains sparse signals. 59 Thus, it is appealing to develop and apply advanced methods that help us to extract useful information from the data for patient classification. The approach proposed integrates genomics data and a melanoma network with an autoencoder model to classify TCGA melanoma patients, and three subgroups with distinct prognostic implications and molecular characteristics were identified. In addition, the complex architectures of deep learning models make them difficult to interpret. 21 We combined the feature-ranking method, the melanoma network and gene set enrichment analysis to interpret the genomic features learned by the autoencoder model for patient stratification. Such a method not only allowed us to provide the biological explanation of the genomic features at the pathway level but also to characterize their molecular features such as mutation and interactome profiles.
We used a melanoma network that is manually constructed and curated using the literature and further expanded using relevant databases. Such a knowledge-based network allowed us to put the newly discovered relevant interactions into the context of the existing knowledge and facilitate the mining and interpretation of the genomics data. However, it may be compromised by not considering other intracellular interactions that are known to be important but not yet found to be relevant for melanoma. An all-encompassing network can be reconstructed using databases such as BioGRID 60 and Omnipath 61 but may introduce redundant information and thus increasing the computational cost and potential risk of overfitting while training the deep learning model. Furthermore, there is no guarantee that an identified community structure of a network is unique because the procedure is complicated by a multiplicity of community detection algorithms and multiple and conflicting definitions of communities. 62 We identified several different partitions with comparable modularity values, and systematic analysis of their properties ensured that the selected community structure of the melanoma network is optimal and representative.
The autoencoder model transforms the genomics data of the identified network communities to a score separately, and one could take advantage of the fact that the melanoma network is a whole and design a model in which all the network community scores are simultaneously learned. However, although such a model may be biologically more reasonable, it may face a trade-off that requires more computation time and increases the number of parameters such as weights and biases, therefore raising the risk of overfitting due to the high dimensionality of genomics data. Furthermore, applying deep learning approaches to high-dimensional genomics data faces the interpretation challenge. Model interpretation is often more important than predicting with high accuracy, for example, identifying prognosis factors for predicting the survival time of patients. 10 However, hidden nodes, computed by fully connected layers, are not able to represent explicit biological components. We demonstrated the interpretable ability of the approach, combining SHAP analysis, annotated pathway knowledge and the melanoma network, to provide contributing features for patient classification with the biological explanation.
In the future, we can modify the model by linking to survival prediction 63 or considering multiple tasks. 64 The upgraded model may serve as a novel, powerful tool for dermato-oncologists to identify patients with a high risk for progressive disease and poor survival who will benefit from modified, cluster-adapted follow-up schemes or more aggressive treatment strategies. Vice versa, patients with a low risk according to their gene expression and DNA methylation profiles may also take advantage of the model by saving them from unnecessary treatment-related toxicities. They may receive "wait-and-see" care only instead of adjuvant treatments that are often accompanied by severe adverse events and decrease the patients' quality of life.
Therefore, such kinds of deep learning model may support physicians in selecting the best individual therapy and follow-up scheme for each patient based on genomics data.
To date, the usability of the model is limited as generating gene expression and DNA methylation data from single patients remain costly and currently not feasible in a daily clinical setting. Additionally, processing, analysis and interpretation of sequencing data are still challenging. However, the model has the potential to facilitate analysis and interpretation by reducing a large amount of genomics data to network community scores for stratifying patients. Furthermore, the model may serve as a powerful tool for complementing already existing predictive biomarker panels that investigate the expression of a few genes, for example, the MelaGenix score. 65 Besides, a future step might be the validation of the model using an independent patient cohort and expanding such a model on treatment response to immunotherapy.
Taken together, the approach provides an effective integration of genomics data and a disease network into a deep learning model.
Although we demonstrated the capability of the approach with only melanoma data, it is a generic approach and therefore applicable to other cancer types and possibly other types of diseases, too. Besides, our study also shows the necessity of developing and implementing multiple analysis approaches given the complexity of multidimensional data.

ACKNOWLEDGMENTS
We thank the colleagues at the Department of Dermatology for their valuable input and discussion. Due to the space constrain of the journal, we, unfortunately, cannot cite all relevant publications, especially the used R packages and experimental evidence for elaborating the impactful features identified by the analyses. Open Access funding enabled and organized by Projekt DEAL.

CONFLICT OF INTEREST
The authors declare no potential conflict of interest.