Survival prediction in breast cancer remains a complex task due to the intrinsic high dimensionality, noise, and heterogeneity of multi-omics datasets, which pose significant challenges for conventional predictive models. Existing machine learning methods often fail to fully leverage the complementary information embedded across different omics layers and struggle to generate clinically interpretable outputs. In this study, we propose a hybrid BiLSTM+CNN deep learning architecture that effectively captures both temporal dependencies and hierarchical feature representations within integrated multi-omics data. The model demonstrates superior predictive performance on benchmark METABRIC and TCGA datasets, while incorporating interpretability mechanisms to enhance clinical relevance. By addressing both the data integration and interpretability bottlenecks, this work provides a robust and scalable framework for precision oncology applications, offering improved survival prediction capabilities that can directly inform personalized treatment strategies.

This study uses the METABRIC breast cancer dataset, consisting of 1980 patient records, available through the cBioPortal database ( https://www.cbioportal.org/study/summary?id=brca_metabric). ¹³ The cBioPortal offers a web-based platform for exploring and visualizing multidimensional cancer genomics data, converting complex molecular profiling from cancer tissues and cell lines into readily understandable genetic, epigenetic, gene expression, and proteomic information. The dataset contains information from three data modalities: clinical profile, gene expression profile, and copy-number alteration (CNA) profile. Patients were grouped based on their survival outcomes into two categories: long-term survivors (≥5 years) with 1489 samples (labeled as ‘0’), and short-term survivors (<5 years) with 491 samples (labeled as ‘1’). The median age at diagnosis for patients is 61 years, with an average survival duration of 125.1 months. Table 1 summarizes the key characteristics of the METABRIC dataset.

The clinical features available in the METABRIC dataset include age at diagnosis, tumor size, estrogen receptor status, HER2 status, and stage at diagnosis. During pre-processing, two of the original 27 clinical features were removed due to missing data and redundancy, reducing the number of clinical features to 25. This feature reduction ensures that only the most relevant variables are retained, enhancing the predictive capacity of the proposed model. To provide deeper insights, a univariate t-test analysis was conducted on key clinical features to assess how they differ between short-term and long-term survivors. The results of the descriptive statistics and t-tests are presented in Table 2 above. Although the METABRIC dataset contained 25 clinical features after preprocessing, the t-test analysis highlighted a subset of widely recognized prognostic indicators, including age, tumor size, estrogen receptor status, HER2 status, and stage at diagnosis. These variables were chosen because they are consistently associated with breast cancer survival outcomes and provide clinically interpretable insights, as also emphasized in prior studies. ¹⁴ The remaining clinical variables were retained for model training but are not individually reported here, as they exhibited limited statistical differentiation between survivor groups.

As shown in Table 2, several clinical features differ significantly between the two survivor groups. Long-term survivors tended to be older at the time of diagnosis, with a mean age of 63.2 years, compared to 55.3 years for short-term survivors (p = 0.03). In addition, tumor sizes were notably smaller among long-term survivors (2.8 cm) compared to those in the short-term group (4.2 cm), with a p-value of 0.01. Similarly, the stage at diagnosis was lower for long-term survivors, indicating early detection and more favorable prognoses (p = 0.02). However, HER2 status showed no statistically significant difference (p = 0.12), suggesting it may not directly influence survival outcomes. These findings emphasize the role of early detection, tumor size, and diagnosis stage in predicting long-term survival, which are essential factors for clinical decision-making.

The heterogeneity in the METABRIC dataset reflects the merging of data from multiple hospitals, leading to variability in clinical practices, treatment protocols, and laboratory standards. In addition, patients often received concomitant medications alongside primary treatments, such as supplements or medications to manage side effects. These factors underscore the complexity of predicting survival outcomes in breast cancer patients. To further validate the model’s generalizability, the TCGA-BRCA dataset was employed. This dataset, available through the GDC portal ( https://portal.gdc.cancer.gov/projects/TCGA-BRCA ), contains 1080 patient records with the same three data modalities as METABRIC—clinical profile, gene expression profile, and copy-number alteration (CNA) profile. ¹⁵ Table 3 below summarizes the details of the TCGA-BRCA dataset.

Deep learning (DL) models have demonstrated remarkable achievements in tasks involving histological images, topographies, and clinical data. However, their performance with gene expression data remains constrained due to the complex nature and high dimensionality of such datasets, which often require thousands of instances to achieve reliable outcomes. To mitigate these challenges, Data Augmentation (DA) techniques can be adapted for transcriptomic data, although their application is less common than in imaging tasks. In this study, we employed noise injection as augmentation strategies. Noise injection has been shown to be particularly effective for omics data, as it simulates the measurement variability inherent in RNA-seq and microarray platforms while generating synthetic yet biologically plausible training instances. Prior studies such as Islam et al. demonstrated that injecting noise into gene expression features improved classification robustness across cancer datasets. ¹⁶ Additionally, gene expression data introduces challenges due to its variability and susceptibility to the curse of dimensionality—where datasets contain more features than available samples. As a result, DA techniques play a pivotal role in increasing the size of training datasets by generating synthetic data samples, thereby improving the generalization capacity of models.

In this study, random rotation and noise injection techniques were applied to gene expression data as DA methods. The noise injection technique involved randomly selecting training samples and altering up to 25% of their features. The noise was generated from a normal distribution with a standard deviation of 0.2 and was added to the original feature values. To ensure data validity, the modified values were clamped within the range of [0, 1]. The selected standard deviation value (0.2) ensured that the augmented samples remained close to the original data instances. Random rotation, on the other hand, was adapted to the feature space rather than physical geometry. Each gene expression profile was treated as a high-dimensional vector, and random orthogonal transformations were applied to rotate these vectors while preserving their variance and overall structure. This produced synthetic samples that retained the statistical properties of the dataset but in alternative orientations of the feature space. Such transformations are conceptually related to PCA-based rotations and latent-space perturbations used in omics DA studies. ^{16 , 17} From a biological perspective, these augmentations capture natural variability in gene–gene correlation structures across patients or cohorts, helping the model learn robust interaction patterns instead of relying on fixed dependencies. A similar idea was demonstrated by previous study, where latent-space mixing in single-cell RNA-seq generated biologically plausible samples that enhanced generalization. ¹⁷ Recent research on DL models for genomic datasets highlights the potential benefits of DA techniques, although the application of DA to genomic data remains relatively unexplored. ^{18, 19} The integration of DA techniques in this study addresses imbalances in gene expression data and enhances the predictive capabilities of the model by preventing overfitting to limited sample sizes.

To ensure reproducibility of the stochastic processes involved in data augmentation, model training, and dataset partitioning, all random operations were controlled using a fixed global random seed of 42 across NumPy, TensorFlow, and Python’s built-in random module. For model evaluation, ten-fold cross-validation was implemented using sklearn.model_selection.StratifiedKFold (n_splits = 10, shuffle = True, random_state = 42) to maintain class distribution across folds. The same data splits were consistently applied across all modalities and compared algorithms to ensure fair evaluation. Within each fold, the training set was further divided into training and validation subsets using an 80/20 split with random_state = 42.

This study utilized three key data modalities: clinical profile, gene expression profile, and copy-number alteration (CNA) profile. Each of these datasets underwent a thorough pre-processing pipeline to ensure data quality and consistency for analysis. The pre-processing steps involved handling missing values, normalization, feature discretization, and feature selection, all of which were necessary to prepare the data for the deep learning algorithm. Figure 1 shows the flowchart of preprocessing data.

To address missing values in the gene expression and CNA datasets, the weighted nearest neighbor (KNN) algorithm was employed. ²⁰ This algorithm identifies the nearest samples based on Euclidean distance and estimates missing values using a weighted average of neighbouring samples. The imputation process used k = 10 nearest neighbors, with Euclidean distance as the similarity metric and inverse distance weighting applied during estimation. After the imputation, the datasets were normalized to maintain a consistent scale for all features. ⁹ Following normalization, the CNA data was discretized into five categories: −2, −1, 0, 1, and 2, representing varying levels of copy number variation. Similarly, the gene expression values were categorized into three classes: -1 (under-expressed genes), 0 (baseline genes), and 1 (over-expressed genes). These discretization steps ensured that the features were more interpretable and ready for machine learning processing.

Given the high-dimensional nature of both the CNA and gene expression datasets, feature selection was essential to reduce dimensionality and improve the model’s generalizability. The Maximum Relevance Minimum Redundancy (MRMR) algorithm was selected for this task due to its ability to identify features that are highly relevant to the target variable while minimizing redundancy among selected features. Alternative methods, such as LASSO regression and Principal Component Analysis (PCA), were considered; however, MRMR was chosen because it provided better interpretability and reduced the risk of overfitting in our experiments. The MRMR algorithm was implemented using the Mutual Information Difference (MID) criterion to rank and select relevant features, where mutual information was estimated using a discrete mutual information approach. The class variable was defined as a binary survival label (0 representing long-term survival and 1 representing short-term survival). Feature selection was performed using an incremental greedy forward selection strategy, where features were sequentially added based on their MRMR ranking.

Among the three modalities used in this study, copy number alteration (CNA) and gene expression data posed significant dimensionality challenges, with each CNA sample containing over 26,000 features and gene expression profiles exceeding 24,000 features. When the number of features exceeds the number of observations, models tend to overfit easily, a problem commonly referred to as the curse of dimensionality. To address this issue, the feature selection process was conducted using a gradational strategy, where subsets of features were generated in increments of 100 and evaluated using the Area Under the Receiver Operating Characteristic Curve (AUC-ROC) as the performance metric. The MRMR algorithm was executed across feature set sizes ranging from 100 to 500, and the subset yielding the highest validation performance was selected. A 5-fold cross-validation strategy was applied during evaluation to ensure robustness and reduce the risk of overfitting. This tuning process identified 400 features for gene expression and 200 features for CNA as optimal, balancing predictive accuracy with model generalizability. The clinical dataset, which contained 27 initial variables, was reduced to 25 key features after preprocessing. These included well-known prognostic indicators such as hormone receptor status, tumor size, menopausal state, lymph node positivity, histological grade, treatment type, and surgical information, all of which have established relevance to breast cancer survival outcomes.

From the CNA profile, the feature count was reduced from 26,298 to 200, and from the gene expression profile, it was narrowed down from 24,368 to 400. For the clinical dataset, the original 27 features were reduced to 25 after removing two features with missing data. The resulting pre-processed dataset, summarized in Table 4 above, served as the input for the proposed deep learning algorithm, enabling more accurate survival predictions.

Convolutional Neural Network (CNN) is a type of feed-forward neural network widely used for tasks involving image processing, natural language processing (NLP), and time series data prediction. ²¹ One of the key advantages of CNN is its local perception mechanism and weight sharing across different layers. This design significantly reduces the number of parameters, thereby improving the model’s efficiency in training and generalization. A typical CNN model is composed of three essential components: the convolution layer, pooling layer, and fully connected layer. The convolution layers extract relevant features from the input data, though the extracted features may have a high dimensionality. To address this, a pooling layer is applied after each convolution layer, which reduces the feature dimensions and computational cost while retaining the most important information.

While CNN has demonstrated exceptional performance in many domains, it has limited capacity to process large-scale, multi-modal data such as genomic and clinical datasets. Recently, researchers have focused on multi-source data integration to enhance the predictive capabilities of deep learning models. These advanced deep learning algorithms that combine multiple data modalities exhibit superior performance over models that rely solely on a single data source. The integration of CNN into such frameworks makes it a promising tool for tasks like cancer survival prediction by efficiently capturing feature patterns across different data types.

The traditional Long Short-Term Memory (LSTM) network, though effective in modeling sequential data, processes information in only one direction—either forward or backward in the sequence. ²² This limitation can hinder the model’s ability to fully capture the temporal dependencies inherent in sequential datasets. To overcome this, the Bi-directional Long Short-Term Memory (BiLSTM) network was developed, enabling the processing of information in both directions—forward and backward. The core idea of BiLSTM is to analyze sequences both front-to-back and back-to-front. In this model, one LSTM layer processes the sequence from the start to the end, while another layer processes it from the end to the start. This dual-directional processing allows the network to retain information from both past and future contexts, making it particularly useful for analyzing time series data and sequential inputs from multi-omics datasets.

In this study, the input data from multi-omics sources is first processed by two BiLSTM layers in the initial module. The extracted features are then passed to the CNN layers in the subsequent module for further feature extraction and dimensionality reduction. The combined BiLSTM and CNN architecture ensures that both the temporal dependencies and spatial patterns in the data are captured. In the final stage, the fully connected layers generate predictions about patient survival, classifying breast cancer patients as either short-term or long-term survivors.

This study presents an improved deep learning algorithm by integrating a Bi-directional Long Short-Term Memory (BiLSTM) network with a Convolutional Neural Network (CNN) to predict breast cancer survival and extract meaningful features from multi-omics data. The BiLSTM addresses the limitations of traditional LSTM networks, while CNN complements it by capturing the spatial patterns in the data. The proposed combination offers superior performance by leveraging the strengths of both models: BiLSTM for temporal sequence learning and CNN for feature extraction. An overview of the proposed model is illustrated in Figure 2 below:

During the initial phase, both BiLSTM and CNN layers are configured with specific filters to extract key features. These features are processed through convolution and dense layers, generating a feature map that feeds into subsequent stages of the model. As described in, ¹⁴ the Glorot normal initializer is used to initialize the filter values, ensuring that the parameters follow a normal distribution with a mean of zero. A fixed seed value of 0.1 is used to maintain consistency in model training, preventing variation in results between different runs.

Hyperparameter tuning was conducted using a grid search approach. The key parameters tuned included the number of layers, filter sizes, learning rate, and regularization strength. For each parameter combination, model performance was evaluated on the validation set using AUC-ROC as the primary metric. The final configuration (outlined in Table 5) was selected based on the highest AUC and accuracy scores observed during cross-validation. This tuning process ensured optimal performance while preventing overfitting. The parameters and architecture of the BiLSTM+CNN algorithm is detailed in Table 5 below:

Below is an overview of the BiLSTM+CNN algorithm:

The pre-processed data serves as input for the BiLSTM layers, whose outputs are subsequently passed to the CNN layers for further feature transformation. Within the convolutional layers, a stride size of 2 is employed, meaning the filter shifts by 2 units across the input matrix during the convolution operation. Padding is applied to ensure that feature sizes remain consistent. A flattened layer follows, converting the multi-dimensional output into a format suitable for the dense layer, which comprises 150 hidden units. To mitigate overfitting, L2 regularization is incorporated within the CNN model. ²³ The activation functions used include ReLU for the convolutional layer and sigmoid for the dense layer. The Adam optimizer ²⁴ was employed for optimization, with binary cross-entropy serving as the loss function for the binary classification task of predicting patient survival.

The proposed BiLSTM+CNN algorithm is divided into three phases: •

Phase 1: The algorithm is trained using clinical data, CNA data, and gene expression data.

Stacked RF model is depicted in Figure 3 below:

The final classification employed a Random Forest (RF) model with n_estimators = 200, max_depth = None, random_state = 0, and class_weight = ‘balanced’. These parameters were selected through grid search, ensuring the best AUC and accuracy scores during cross-validation. This configuration offered the optimal balance between precision and recall, especially for the imbalanced classes in our dataset. Using class_weight = ‘balanced’ mitigated the risk of overlooking minority classes, while setting max_depth = None enabled the model to capture complex feature interactions without overfitting.

The performance of the proposed BiLSTM+CNN algorithm was evaluated using several metrics, including Sensitivity, Specificity, Precision, and Accuracy and the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) curve. The metrics are defined as follows: Sensitivity = Tp Tp + Fn (1) Specificity = Tn Tn + Fn (2) Precision = Tp Tp + Fp (3) Accuracy = Tp + Tn Tp + Tn + Fp + Fn (4)

Here, TP (true positive), TN (true negative), FP (false positive), and FN (false negative) denote the classification outcomes. Additionally, the AUC-ROC curve assesses the model’s ability to distinguish between classes across various thresholds, providing a comprehensive view of performance beyond a single point metric.

The ten-fold cross-validation approach was adopted for model evaluation, following recommendations from prior studies. ¹⁴ In this method, the dataset is randomly divided into ten equal subsets. For each fold, nine subsets are used for training, while one subset is held out for testing. This process ensures that every data point is used for both training and testing, thereby providing a more reliable performance estimate.

Within each merged training set, 80% of the data is allocated for training the model, while the remaining 20% is reserved for validation to fine-tune hyperparameters and prevent overfitting. The Keras and TensorFlow libraries were employed for model implementation, ensuring computational efficiency and ease of experimentation.

The computational environment, all experiments were conducted on a workstation running Ubuntu 20.04 LTS equipped with an NVIDIA GeForce RTX 3090 GPU (24 GB VRAM) and 64 GB RAM. The software stack comprised Python 3.9.7, TensorFlow 2.8.0, Keras 2.8.0, scikit-learn 1.0.2, NumPy 1.21.5, and pandas 1.3.5. A configuration file and step-by-step reproducibility protocol are provided in the archived software repository ( https://doi.org/10.5281/zenodo.15964646).

The proposed deep learning algorithm leverages BiLSTM and CNN for feature extraction from multi-omics data. The AUC metric from the ROC curve, along with accuracy, is used to evaluate the model’s performance. Figure 4 below shows the ROC curves of the BiLSTM+CNN compared to CNN for the METABRIC dataset. The AUC values are 0.90, 0.87, and 0.87 for clinical data, CNA, and gene expression data, respectively. To provide a more comprehensive view of the model’s performance, we report the 95% confidence intervals (CI) for each modality. These CIs offer an estimate of variability across different trials, ensuring more reliable interpretation of results: •

Clinical data: Accuracy = 0.90, 95% CI [0.9027, 0.8973]

These results indicate the model’s robustness across different data modalities, though performance on gene expression data is slightly lower, reflecting the challenges posed by high-dimensional data. The third plot in the above Figure 4 is presented with data points rather than a continuous line to highlight specific thresholds along the ROC curve. This visualization helps illustrate how key decision points, such as cutoff thresholds, impact the true positive (TPR) and false positive rates (FPR). Disconnected lines may occur due to discrete prediction values or gaps in the input data, especially when thresholds do not span the full range of possible values. While this approach enhances interpretability, future iterations could explore smoothing techniques to provide a continuous curve.

Our model’s false positive rate was higher than expected, which could have clinical implications. To address potential overfitting and variance due to the small dataset size, we applied ten-fold cross-validation. The 1980-patient dataset was split into ten subsets, with nine subsets for training and one for testing. Each training set was further divided into 80% for training and 20% for validation.

We combined extracted features from BiLSTM and CNN into a stacked feature set, which was then classified using a Random Forest (RF) algorithm. As previous studies show, RF performs better with stacked features compared to other classifiers. ¹⁰ Performance metrics, including sensitivity, specificity, and precision, were calculated to assess the model’s effectiveness.

Figure 4 above presents the ROC curve comparing BiLSTM+CNN and CNN feature extractors on METABRIC data. The results demonstrate superior AUC values for BiLSTM+CNN across different modalities. Table 6 below summarizes the comparison of AUC and accuracy with existing algorithms:

The bold values in Table 6 highlight the best performing results for each data modality (clinical, CN and gene expression) across all compared algorithms. The results clearly demonstrate that the proposed BiLSTM+CNN algorithm performs better than previous algorithms across multiple data modalities. However, in the gene expression modality indicate that the SiGaAtCNN Stacked RF model achieved superior performance, outperforming the proposed BiLSTM+CNN model in that specific category. The comparison included models such as MDNNMD, SiGaAtCNN, and Heterogeneous Stacked RF. As shown in Table 7 above, our algorithm outperforms others in terms of accuracy, precision, sensitivity, and Matthews correlation coefficient (MCC).

To further validate the performance, we used the TCGA-BRCA dataset. ¹⁵ This dataset contains 250 long-term survivors and 830 short-term survivors, with data modalities matching those in the METABRIC dataset. Pre-processing was conducted using the same steps outlined in Sections B and C.

Figure 5 presents the ROC curve for the TCGA dataset, demonstrating that the BiLSTM+CNN Stacked RF algorithm maintains high performance across datasets. Below are the performance metrics along with the 95% confidence intervals (CI): •

Clinical data: Accuracy = 0.739, 95% CI [0.741, 0.737]

These results show that the model generalizes well to the TCGA dataset, especially on gene expression data, where it achieves high accuracy. The 95% CI for each modality further supports the robustness of the proposed model. Despite inherent differences between the METABRIC and TCGA datasets, the BiLSTM+CNN model achieves high accuracy across all data modalities. The results of our algorithm are compared with other state-of-the-art algorithms in Table 8 as follows:

Our study builds upon existing research that utilizes multi-omics data for survival prediction. Curtis et al. ¹³ identified prognostic biomarkers using a multidimensional competition-based framework with the METABRIC dataset, while our study advances this work by integrating BiLSTM and CNN architectures for capturing temporal and spatial patterns across data modalities. Unlike the framework by Curtis et al., ¹³ which focused primarily on identifying subgroups, our model emphasizes multi-omics data integration for improved survival predictions and interpretability.

Additionally, Yousefi et al. ⁷ and Mobadersany et al. ¹⁸ employed convolutional networks for cancer survival outcome predictions, but their models primarily focused on histological data. In contrast, our model integrates clinical, CNA, and gene expression data, providing a more comprehensive and interpretable prediction framework. This integration allows the model to extract complex patterns that go beyond histological data alone. The study by Jadoon et al. ²⁵ proposed a heterogeneous multiple kernel learning approach for breast cancer prognosis, addressing the challenge of multimodal data. While their approach is robust, our deep learning-based solution offers enhanced predictive performance through the combined use of BiLSTM and CNN architectures, which capture both sequential and spatial information across data types. Similarly, Phan et al. ⁶ demonstrated the use of machine learning models for decoding breast cancer with multi-omics data but faced challenges related to model interpretability and high dimensionality. Our approach, with MRMR feature selection, addresses these challenges by reducing dimensionality while retaining the most informative features.

In summary, our model offers a novel combination of deep learning models and feature selection techniques to provide actionable clinical insights. The use of decision-level integration ensures robust predictions across datasets, with significant improvements observed on both METABRIC and TCGA datasets. These comparisons highlight how our work advances the field by building on previous methodologies while addressing key limitations, such as the interpretability and scalability of predictive models.