Predicting Melanoma Immunotherapy Efficacy: Neural Network Models with Gene Expression and Clinical Data

1. Introduction

Recent research has made significant advancements in treating melanoma through Immune Checkpoint Blockade (ICB). Within cancer cells, checkpoints regulate the strength of the immune response, preventing excessive attack by T cells. Blocking these checkpoints with specific antibodies stimulates the immune system's response against cancer cells. Therapies targeting checkpoint molecules such as Cytotoxic T Lymphocyte Antigen 4 (CTLA4), Programmed Cell Death 1 (PD1), and PD1 ligand 1 (PDL1) have shown early success in clinical trials. Approved treatments now include anti-PD1 drugs like nivolumab or pembrolizumab, and anti-CTLA4 therapy like ipilimumab, either alone or in combination. However, not all patients benefit from these treatments, prompting extensive research into the molecular mechanisms of resistance to ICB. In addition to tumor-intrinsic and tumor-extrinsic biomarkers, patient genomic characteristics are also considered when making treatment decisions. Studies based on Ribonucleic Acid (RNA) analysis have identified gene expression patterns associated with immune cell infiltration within the tumor microenvironment, which in turn correlates with patients' overall survival. Specifically, expressions of Human Leukocyte Antigens (HLA)-A and Cluster of Differentiation (CD)8A show promise as potential targets for CTLA4-directed therapies. Some studies suggest that neoantigens and viral presence can reveal mutations enabling tumors to resist immune attacks, although their utility as predictive markers remains uncertain. Furthermore, prior failure of Mitogen-Activated Protein Kinase (MAPK)-targeted therapy negatively impacts subsequent responses to ICB in melanoma patients.Top of Form

Recent findings implicate that Deoxyribonucleic Acid (DNA)-and RNA-level genomic information has predictive value for the clinical benefit of both anti-CTLA4 and anti-PD1 therapy. For example, while having many genetic mutations might indicate longer overall survival, it doesn't always mean a better response to anti-PD1 therapy. Therefore, there's a need to explore other genetic and non-genetic factors influencing treatment outcomes. To address this, extensive analysis of large patient groups is necessary. Finding ways to predict treatment response has been challenging due to the complexity of the immune system's response, differences in tumor environments, and the vast amount of data generated. To handle this data, sophisticated and automated techniques are required. Researchers have developed statistical models, like those by [1], which consider associations between genetic and transcriptomic features and treatment response. Because transcriptomic data often contains many interrelated variables, these models must account for nonlinear relationships. Recent studies have shown that Deep Neural Networks (DNNs) [2,3,4] can effectively model these nonlinear relationships, leading to better predictions of cancer therapy response. Additionally, neural network-based approaches have demonstrated superior performance in classifying gene expression samples such as [5].

Despite the growth of Deep Learning (DL) methods, their application in predicting immunotherapy response remains limited. However, there's a growing need for such methods, especially in melanoma patients undergoing ICB therapy. Artificial Neural Networks (ANNs) serve as the backbone of these predictive models. Using next-generation sequencing data, we aim to predict individual responses to drug therapy. However, the high-dimensional nature of the data and the limited number of samples present challenges. Therefore, we follow established scientific approaches, including quantifying drug response, selecting molecular features, fitting DL models, and evaluating their performance. In our study, we focus on quantifying Progression-Free Survival (PFS) and Overall Survival (OS) as measures of drug response. We also consider other genomic and clinical features, such as tumor mutational burden and patient demographics. Moreover, we account for prior treatments, like MAPK inhibition, if available. By developing and evaluating DL models, we aim to provide clinicians with sophisticated tools for treatment decision-making. However, working with high-dimensional data poses several challenges. Overfitting, where models perform well on training data but poorly on new data, is a common issue [6,7,8,9,10,11,12]. To address this, we employ various techniques such as feature selection, regularization, and data augmentation. Additionally, the imbalance in class distribution within our datasets requires special attention, which we address using methods like Synthetic Minority Over-sampling TEchnique (SMOTE) [13,14,15]. Validating our results is crucial but hindered by limited data availability. Therefore, we emphasize the importance of data sharing and standardization. Despite these challenges, integrating data from multiple sources can enhance the predictive power of our models. While our focus is on PFS/OS, our approach can be adapted for other endpoints or treatment modalities. The objective of this study is to automate the prediction of immune checkpoint therapy success in melanoma patients using gene expression data alongside genomic and clinical features, with patients as the primary outcome variables, addressing high-dimensional data challenges and small sample sizes inherent in DL.

The study is as follows; similar papers are shown in the following section. The materials and methods are provided in Section III. The experimental analysis is carried out in Section IV. The discussion is given in Section V, after which we draw some conclusions and outline our plans for further research in Section VI.

2. Related Works

Recent advances in the treatment of melanoma have demonstrated significant advances in ICB. Laboratories in cancer cells regulate the immune system, which regulates the action of T cells against cancer cells. Monoclonal antibodies that interfere with these checkpoints have shown promise in the immune response against cancer cells. Remarkable success has been observed with antibodies targeting checkpoint molecules such as [16,17]. Approved therapies include anti-PD1 (nivolumab or pembrolizumab) [18] and anti-CTLA4 (ipilimumab) [16], individually or in combination [19]. However, the efficacy of these treatments varies between patients, prompting extensive research into the molecular mechanisms of ICB resistance [1]. In addition to tumor intrinsic and extrinsic biomarkers, genomic profiling of patients is being investigated for treatment stratification such as [20]. RNA-based studies have identified gene expression signatures of immune penetration within the tumor microenvironment, which correlate with patient survival such as [21,22]. Whereas, this study highlights the expression of their potential CTLA4-targeted antibodies, such as HLA-A and CD8A. However, challenges remain in the search for predictive responses due to the complexity of the immune system and heterogeneous tumor states. Addressing these challenges requires sophisticated research approaches, with predictive models more flexible [1] and deeper tissues [23]. In response prediction modeling complex relationships have shown promise, but their role in ICB treatment of melanoma has been largely unexplored [5]. Given the demand for immunotherapy response prediction methods, especially for melanoma patients undergoing ICB therapy, the integration of DNNs-based models could significantly improve treatment decisions such as [23]. He emphasizes that, such models provide an effective way to capture the complexity of cancer biology compared to traditional linear models. Therefore, the search for effective prediction models in predicting immune responses, especially the use of DL techniques, presents an important frontier in melanoma research, with great potential for personalized treatment approaches.

3. Materials and Methods

3.1. Data Analysis

In this study, we utilized Boston property price 1 dataset which is encompassed within the Sklearn package. It includes 506 instances with 13 numerical/categorical attributes and the target variable. These attributes comprise details such as locality, surroundings, and crime incidence. The dataset has been utilized to tackle a regression issue by forecasting the median worth of residences. The original data is accessible in the Sklearn manual. All gene expression profiles employed in this investigation were obtained from published melanoma assemblies treated with diverse forms of ICB therapies. Due to differences in size and regimen among datasets, it is imperative to categorize the assemblies based on therapy type and underscore the size of each subset. Within each assembly, a patient may have received either one immune checkpoint antagonist (monotherapy), one antagonist consecutively (sequential therapy), or a blend of antagonists (combination therapy). Furthermore, certain melanoma patients have previously undergone MAPKi therapy. For each dataset, the Transcripts Per Million (TPM) values of gene expression along with data regarding treatment and survival time were extracted from the published collections. Depending on the assembly, supplementary details about patients or tumor milieu were also accessible. Nonetheless, such supplementary information is frequently limited to a single assembly, complicating the amalgamation of assemblies. For instance, age is documented in the [24]—Gide/[18]—Hugo/[25]—Van Allen assembly but not in the [1]—Liu/[26]—Riaz assembly. Given that the Liu dataset, comprising 121 records, is the most extensive and the efficacy of DL techniques aligns with the acquisition of sizable, standardized datasets, this dataset was predominantly employed to assess the ANN methodology. The Liu dataset was procured from the [1]. Patients within this assembly were administered either pembrolizumab or nivolumab as monotherapy or ipilimumab and pembrolizumab/nivolumab as sequential therapy. Some patients also received antecedent MAPKi intervention. Alongside the 20,849 gene expression values in TPM format, six clinical attributes were selected for analysis: gender, heterogeneity, nonsynonymous mutation, cancer stage, elevated Lactate Dehydrogenase (LDH), and treatment type (pembrolizumab or nivolumab). The Gide dataset was sourced from the [24]. The individuals within this dataset have received treatment either through single-agent therapy employing pembrolizumab or nivolumab, or through combined therapy by incorporating ipilimumab as the primary agent, respectively. Patients harboring BRAFV600 mutations were classified as individuals with previous MAPKi exposure. The Riaz dataset was acquired from the [26]. This particular group comprises patients subjected solely to nivolumab therapy. There was no available data regarding prior MAPKi treatment. The Hugo dataset was obtained from the [18]. This set comprises patients treated with pembrolizumab monotherapy. Patients with previous MAPKi exposure were included in this cohort. The Van Allen dataset was sourced from the [25]. This dataset involves patients undergoing ipilimumab monotherapy. Patients who received a BRAF inhibitor subsequent to ipilimumab treatment were categorized as individuals with prior MAPKi exposure.

3.2. Model Analysis

The described ANN process is documented using the publicly accessible web application Jupyter Notebook in conjunction with the Python programming language. Numerous Python packages were utilized for constructing the ANN workflow. These packages were installed via pip (20.3.3) and deployed on Debian GNU/Linux 10 (buster). The workflow can be segmented into four primary phases. Initially, the data necessitates preprocessing and division into training and testing datasets. An optional stage follows, concerning dimensionality reduction through feature selection. Subsequently, the ANN is trained and tested, encompassing parameter tuning with GridSearch, ANN training with optimized parameters, and final model validation using the held-out test dataset. Data preprocessing emerges as a pivotal aspect of the ANN modeling workflow. Initially, all samples featuring missing values in any attribute were eliminated, resulting in 506 samples for the Boston Housing dataset and 118 for the Liu dataset. Given that the dataset comprises both numerical and categorical variables such as treatment, all input variables necessitate conversion into numerical representations. In the Boston Housing dataset, all 13 attributes are inherently numeric, while categorical attributes in Liu were transformed utilizing the LabelEncoder function from Sklearn. For binary predictions, outcomes were stratified using an appropriate threshold value. In the Boston housing dataset, a MEDV (Median value of owner-occupied homes in $1000) less than 22 is categorized as low-price houses, and greater than 22 as high-price houses. In Liu, patients with PFS of less than 365 days or OS of less than 730 days are classified as short survivors, while those with PFS exceeding 365 or OS surpassing 730 are considered long survivors.

Addressing imbalanced data, the SMOTE-NC oversampling technique is applied to divide the data into two groups, generating synthetic data using a k-nearest neighbor algorithm [27,28,29,30,31,32]. Feature scaling is conducted to standardize feature scales and enhance the gradient descent algorithm's efficiency. Standardization, facilitated by the StandardScaler method from Sklearn, is preferred over Min-Max-Scaling due to its resilience to outliers. Data scaling is solely computed on the training set to prevent data leakage. Recursive Feature Selection (RFE) from Sklearn is employed for feature selection, utilizing Support Vector Regression (SVR) for continuous prediction or Support Vector Classification (SVC) for binary classification. ANN architecture comprises an input layer accommodating gene expression TPM values and clinical variables, followed by one or more hidden layers. Weights and biases initialization employs 'glorot_uniform' and 'zeros' respectively, while training employs the Adam optimizer with variable learning rates. The hidden layer incorporates dense layers followed by dropout layers, utilizing ReLU activation functions for input layers, and variable L₁ and L₂ regularization for kernel regularization. The number of neurons and hidden layers are adjustable parameters. The output layer of the neural network comprises a single node, with activation and loss functions tailored based on the nature of predictions—sigmoid activation and binary_crossentropy loss for binary outcomes (e.g., predicting PFS over or under 365 days), and linear activation with Mean Squared Error (MSE) loss for continuous outcomes (e.g., predicting exact PFS days). The optimizer employed is Adam, with additional variants such as Stochastic Gradient Descent (SGD) with momentum and Nesterov momentum. Several parameters are fine-tuned to optimize the neural network's performance. These include the number of neurons per layer, the number of hidden layers, the learning rate, the L₁ regularization rate, the L₂ regularization rate, the type of optimizer, and the dropout rate. The data is initially split into 80% for training the model and 20% for testing the model's performance on unseen data. Only the 80% training data (n = 94) are used for tuning and training the model. Validation is performed with the remaining 20% (n = 24), as described earlier. The neural network's variable parameters undergo fine-tuning via GridSearchCV with 5-fold cross-validation and a 20% validation split. Once optimal parameters are identified, training commences in shuffle mode to ensure data independence. Two Keras callbacks, EarlyStopping and ReduceLROnPlateau, expedite training and optimize the learning rate. EarlyStopping ceases training if validation loss stagnates for 20 epochs, while ReduceLROnPlateau adjusts the learning rate by a factor of 0.1 if validation loss stalls for five epochs until reaching a minimum of 1e^-6. To evaluate the trained ANN, multiple 80/20 splits are executed using a ShuffleSplit cross-validator, employing tailored metrics. For binary classification, these metrics encompass sensitivity, specificity, positive and negative predictive values, and the confusion matrix. Conversely, for continuous outcomes, evaluation metrics include Mean Absolute Error (MAE), Root Mean Squared Error (RMSE), and R².

4. Experimental Analysis

4.1. Sanity Check with Boston Housing Dataset

To assess the predictive ability of the neurons under ideal conditions, we performed experiments using the well-known Boston Housing dataset as shown in Table 1. This list is typically used to predict housing prices based on 13 distinct characteristics. To simulate the binary prediction task, we divided the target variables into two groups: "affordable" and "affordable" building prices, using a threshold value of 22 (equivalent to $22,000). Parameters chosen for GridSearch optimization include: neuron count = 8,16,32,64 respectively; SecretPopulation = 1,2,3; Best performer = SGD, Adam; study rate = 0.01; Release = 0.0,0.2,0.4,0.6; L1 Constants = 0.0,0.2,0.4; L2 constant processing rate = 0.0,0.2,0.4. The constructed workflow achieved up to 95% accuracy for five train/test splits typically 80/20. However, the prediction performance of the workflow varied depending on the choice of fragmentation, resulting in false predictions and false negatives, and in addition, the specific train/test fractionation had the best coefficient affected by the grid search method. Notably, the number of hidden neurons ranged from 1 to 3 and the number of neurons per layer ranged from 32 to 64 neurons. For a more detailed analysis and considering that the Boston housing dataset is designed for regression functions rather than binary classification, we evaluated the predictive performance of the project using continuous effects. For five different train/trial separations, the root network produced a RMSE of 3.41, MAE of 2.94, R-class, with the same GridSearch parameters used. Furthermore, the optimal parameters differed depending on the specific train/test split, although no constants were used for any of the five splits.

4.2. Predictions on Liu Dataset

The current section present the application of the developed workflow to address the research query regarding predicting individual responses to therapy based on gene expression, in conjunction with described genomic and clinical features. The overarching objective of guiding therapy decisions through a patient's genetic makeup necessitates rigorous validation of predictive models concerning the outcome of ICB therapy based on individual genomic profiles. Given the pivotal role of extensive and standardized datasets in DL success, initial preference is given to the Liu dataset due to its extensive RNA sequencing data as shown in Table 2. This dataset comprises 20849 gene expression values along with six clinical variables, including gender, heterogeneity, nonsynonymous mutation, cancer stage, LDH elevation, and treatment (pembrolizumab or nivolumab). Binary and continuous predictions for PFS and OS were conducted through various splits using identical grid search parameters. The selection of optimal parameters involved a meticulous process, combining testing and knowledge from neural network prediction methods in cancer modeling. Nonetheless, comprehensive details on network configurations, overfitting mitigation strategies, and specific learning parameters were often lacking in available literature. For continuous outcome prediction, addressing the challenge of "exploding gradients" in ANNs proved imperative. This phenomenon arises when significant weight updates during training lead to numerical overflow or underflow. While implementing the gradient norm clipping function in Keras resolved the issue, prediction metrics remained unsatisfactory. Notably, the R² values across all splits were consistently negative, indicating arbitrary deterioration in model performance. Given the complexity of regression tasks for ANNs and associated structural disparities, focus shifted towards binary classification. Moreover, delineating between "long" and "short" survivors proved more clinically relevant in treatment decision-making scenarios. In binary outcome prediction, the model achieved up to 72% accuracy for OS and up to 71% for PFS, accurately classifying patients into respective categories. However, the evaluation of additional metrics such as specificity and sensitivity is crucial, especially considering the modest size of the test split. Upon thorough evaluation, the binary classification model demonstrated promising accuracy rates, with specificity notably outperforming sensitivity in predicting short-term survivors.

To address the challenge of dealing with large amounts of data, a method called feature selection is employed to narrow down the variables to 50 before the actual training of the ANN as shown in Table 3. In the context of predicting PFS, a greater number of features were initially considered. Notably, the selected features demonstrated superior performance in predicting PFS, where none of the initially considered features were present across all five datasets. The feature selection outcomes, detailed in Table 7, highlight consistent identification of genes like HIKESHI and LIN28A across all five datasets, with several others recurring in four out of five datasets. Post-selection, separate training of the ANN for each dataset using these features and clinical variables yielded binary prediction accuracies for PFS ranging from 68% to 86%, surpassing those for OS (49% to 79%). This improvement, alongside a notable boost in PFS prediction accuracy compared to the model without feature selection, underscores the efficacy and adaptability of the approach across diverse dataset characteristics.

4.3. Predictions on Combination of Cohorts

In dealing with the challenge of the limited dataset, we chose to expand the training dataset. To achieve this, we combine the heterogeneous sets, resulting in a pooled dataset of 322 observations. To ensure consistency, we first followed [1], as the TPM cannot accommodate variations in libraries or different levels of rRNA contamination in different clusters. Next, we apply normalization between samples using the mean method of ratios using DESeq2 to reduce systematic and technical biases resulting from different experimental conditions and sequencing protocols. Given the absence of clinical or genomic characteristics beyond OS in all cohorts, our predictions are based solely on gene expression data. Clinical data for each patient in each group are combined using binary coding. For the binary prediction of OS, we train and evaluate in five different partitions using grid search parameters as follows: number of neurons (3, 5, 10), number of hidden layers (1, 2, 3), Number of L₁ regularization (0.0, 0.2, 10). 0.4), the number of L₂ routines (0.0, 0.2, 0.4), and the number of classes (0.0, 0.2, 0.4), including the ADAM optimizer. The predictive accuracy results for OS ranged from 61% to 66% (refer to Table 4). Specificity is consistently greater than 70% for all separations, while sensitivity is high at 55%, indicating better performance in predicting shorter lifetimes than longer ones. Although the metric changes are relatively small compared to the forecasts based on the Liu group alone, overall the combination of the different groups does not significantly improve the forecast results.

5. Discussion

The quest for personalized therapy response prediction poses significant biological and computational challenges, constrained by data accessibility and algorithmic limitations. Utilizing neural networks offers a promising avenue for forecasting therapy responses, particularly based on microarray data. Ensuring the reproducibility of predictive models is paramount for developing clinically relevant techniques, necessitating a flexible neural network architecture and robust data workflow. Applied to the context of predicting responses to ICB therapy, neural networks hold promise for personalized treatment of melanoma patients. However, challenges arise from limited sample sizes in standardized datasets, necessitating innovative approaches to address cohort deficiencies. Notably, diverse treatment modalities further complicate the analysis, highlighting the need for robust computational workflows tailored to biological nuances. Predicting PFS and OS for patients undergoing ICB therapy entails complex analyses integrating gene expression and clinical variables. The developed workflow demonstrates efficacy in leveraging these data, achieving competitive prediction accuracies for both PFS and OS. However, challenges such as overfitting and model variance underscore the importance of regularization techniques and comprehensive parameter optimization. Feature selection plays a crucial role in enhancing prediction accuracy and interpreting model outcomes. Various methods, including RFE and Principal Component Analysis (PCA), offer avenues to address high-dimensional data challenges. Notably, the efficacy of feature selection varies depending on the predictive task, highlighting the need for tailored approaches. Presently, the amalgamation of ICB treatment combinations and limited available data pose challenges to attaining conclusive outcomes, necessitating further investigation, particularly stratifying for treatment combinations. Strides toward enhancing the reproducibility of ICB treatment response prediction have been realized through the formulation of a flexible workflow, the amalgamation of diverse melanoma cohorts predicated on disparate treatments, and the amalgamation of available clinical and genomic features. The adaptability of the workflow renders it amenable to diverse and heterogeneously structured datasets, while the GridSearch function alleviates the human workload by streamlining parameter optimization, fostering clinical applicability. In essence, this study delineates a workflow with demonstrated efficacy provided the available data furnishes pertinent information for regression or classification tasks. By amalgamating several methodologies such as regularization and feature selection, the workflow grapples with high-dimensional, low-sample-size data, which is pivotal for therapy response prediction predicated on transcriptomic data. Moreover, this master study elucidates the challenges stemming from class imbalance and data heterogeneity engendered by therapy prediction grounded on transcriptomic data, underscoring the exigency for further research addressing low-sample-size, high-dimensional data. The evolved workflow, spearheaded in this study, fills the lacuna of DL-based prediction of clinical benefit from ICB therapy for melanoma patients, thus serving as a springboard for personalized treatment decision-making for melanoma patients undergoing ICB therapy.

6. Conclusion and Future Works

Issues persist in predicting response to ICB treatment, largely due to limited melanoma patient samples. Nonetheless, incorporating publicly accessible datasets like those from the Gene Expression Omnibus (GEO) database could bolster predictive models by expanding training sets. Another strategy involves pretraining or transfer learning, where algorithms leverage knowledge from related tasks or unlabeled data to improve performance. Convolutional Neural Networks (CNNs) excel in transfer learning, particularly when faced with challenges in accessing extensive omics data. While this study focused on predicting PFS and OS, exploring alternative endpoints like RECIST criteria presents intriguing possibilities. However, challenges persist in standardizing response measurement for preclinical models, hindering fair model comparisons. Additionally, selecting input features requires careful consideration, with methods incorporating biological knowledge yielding superior results. Enhancing personalized treatment response prediction entails incorporating diverse methodologies like transfer learning and graph-convolutional networks into automated workflows. With access to more data and related datasets, DL has the potential to revolutionize clinical decision-making for ICB treatment in melanoma patients.