Multi-Omic Integration and Machine Learning Reveal Regulatory Networks Driving Breast Cancer Progression

Introduction

Breast cancer remains the most frequently diagnosed cancer and leading cause of cancer mortality among women globally, with an estimated 2.3 million new cases and 685,000 deaths in 2020 (Sung et al., 2021). While significant progress in early detection and targeted therapies has improved survival rates in high-income countries, disparities persist globally, particularly in regions with limited access to molecular diagnostics and personalized treatment approaches (Ginsburg et al., 2020). Traditional breast cancer staging relies primarily on anatomic parameters including tumor size, lymph node involvement, and distant metastasis through the TNM classification system (Giuliano et al., 2017). However, this approach inadequately captures the biological heterogeneity underlying cancer progression, as tumors with identical anatomic staging frequently exhibit markedly different clinical trajectories and therapeutic responses (Yersal and Barutca, 2014).

The molecular revolution in cancer biology has revealed that breast cancer progression from early localized disease to advanced metastatic stages involves coordinated dysregulation across multiple biological layers, including genomic alterations, epigenetic modifications, transcriptional reprogramming, and post-transcriptional regulation (Cancer Genome Atlas Network, 2012). MicroRNAs (miRNAs), a class of small non-coding RNAs approximately 22 nucleotides in length, have emerged as critical regulators of gene expression through messenger RNA (mRNA) degradation and translational repression (Bartel, 2004). Extensive research has documented that miRNAs function as both oncogenic drivers (oncomiRs) and tumor suppressors in breast cancer, with distinct expression patterns associated with tumor initiation, progression, metastasis, and therapeutic resistance (Iorio et al., 2005; O'Bryan et al., 2017; Pencheva and Tavazoie, 2013).

Among the most well-characterized oncomiRs, miR-21 demonstrates pervasive overexpression across breast cancer subtypes and promotes proliferation, invasion, and metastasis through targeting of PTEN, PDCD4, and other tumor suppressors in the PI3K/AKT/mTOR pathway (Bhardwaj et al., 2018). Similarly, miR-155 enhances tumor cell survival and metabolic reprogramming while miR-10b specifically promotes metastatic dissemination (Shi et al., 2015; Ma et al., 2007). Conversely, tumor suppressor miRNAs including the miR-200 family, miR-34 family, and let-7 family undergo frequent downregulation during progression, relieving repression of pro-metastatic factors such as ZEB1, ZEB2, and RAS family members (Gregory et al., 2008; Hermeking, 2010; Johnson et al., 2005). The miR-200/ZEB1/E-cadherin regulatory axis in particular forms a double-negative feedback loop that governs epithelial-mesenchymal transition, a fundamental process enabling cancer cell invasion and metastasis (Burk et al., 2008).

Despite substantial knowledge regarding individual miRNAs in breast cancer biology, their integration with mRNA expression data for clinical classification remains incompletely explored. Previous efforts to develop miRNA-based staging classifiers have reported moderate accuracy (80-85%) using TCGA-BRCA data, substantially lower than the 95-99% accuracy achieved for tumor versus normal tissue classification (Yerukala Sathipati and Ho, 2018; Das et al., 2025). This performance gap likely reflects the non-linear, multifactorial nature of cancer progression, where anatomic stage represents an imperfect proxy for underlying biological state. Recent multi-omic integration studies combining different molecular layers through machine learning have demonstrated 5-10% improvements in classification performance compared to single-platform approaches, suggesting that coordinated analysis of miRNA-mRNA regulatory networks may better capture progression-associated molecular changes (Rappoport and Shamir, 2018; Chaudhary et al., 2018; Qiao and Jin, 2021).

Machine learning approaches offer distinct advantages for multi-omic data integration by identifying complex non-linear patterns and higher-order interactions that traditional statistical methods may overlook (Kourou et al., 2015). Gradient boosting algorithms such as XGBoost have consistently demonstrated superior performance across diverse classification tasks in cancer genomics due to their ability to model feature interactions, handle high-dimensional data efficiently, and provide built-in mechanisms for addressing class imbalance through sample weighting (Chen and Guestrin, 2016; Zhang et al., 2019). Furthermore, ensemble tree-based methods inherently perform feature selection during model training and generate importance scores that facilitate biological interpretation (Sarica et al., 2017).

From a global health perspective, developing robust molecular classifiers for breast cancer progression holds significant implications for clinical decision-making across diverse resource settings. While high-income countries increasingly incorporate multi-gene expression assays such as Oncotype DX and MammaPrint into routine clinical practice, these commercial tests remain cost-prohibitive and logistically challenging in low- and middle-income countries where 60% of breast cancer deaths occur (Yip and Rhodes, 2014; Toss and Cristofanilli, 2015). Identifying minimal informative feature sets from publicly available data that capture core biological progression signatures could enable development of more accessible and affordable diagnostic approaches adaptable to varied healthcare infrastructures.

In this study, miRNA-Seq and RNA-Seq expression data from 1,081 primary breast cancer samples in The Cancer Genome Atlas were integrated to develop and validate machine learning classifiers distinguishing early-stage from late-stage disease. Through comprehensive feature selection, nested cross-validation with rigorous performance assessment, and detailed biological characterization of predictive features, this work aimed to (1) quantify the added value of miRNA-mRNA integration over single-platform approaches, (2) identify and functionally characterize multi-omic biomarker panels associated with progression, (3) establish methodological best practices for addressing class imbalance and optimistic bias in clinical prediction models, and (4) provide biological insights into miRNA-mRNA regulatory networks driving breast cancer progression with potential translational relevance across global populations.

Methods

Data Acquisition and Preprocessing

Clinical annotation, miRNA-Seq, and RNA-Seq gene-level count data for breast invasive carcinoma (BRCA) were retrieved from The Cancer Genome Atlas (TCGA) database through the Genomic Data Commons portal using the TCGAbiolinks R package (version 2.24.0; Colaprico et al., 2016). Sample selection criteria restricted analysis to primary tumor samples with complete clinical staging information and paired miRNA-mRNA expression measurements. From an initial cohort of 1,222 samples, 28 samples lacking complete staging annotations, 67 samples with missing miRNA expression data, 33 recurrent tumors, and 13 metastatic samples were excluded, yielding a final cohort of 1,081 primary tumors (Figure 1a). Samples were stratified into early-stage (Stage I and Stage II, n=822, 76.0%) and late-stage (Stage III and Stage IV, n=259, 24.0%) groups based on pathological stage grouping derived from AJCC 7th edition TNM classifications (Figure 1b).

Sample identifiers were harmonized across data modalities to enable molecular-clinical data integration. Raw miRNA-Seq read count column names contained 'read_count_' prefixes while RNA-Seq columns used TCGA barcodes with different sequencing platform suffixes. Clinical data employed abbreviated 15-character barcodes truncating sequencing batch information. Identifiers were standardized by removing prefixes and truncating to 15-character sample-level root barcodes (participant-sample-vial-portion-analyte), enabling matching of 1,081 samples with complete tri-modal information. The complete analysis workflow from data acquisition through network analysis is illustrated in Figure 1c.

For expression normalization, standard bioinformatics practices adapted for count-based sequencing data were applied (Robinson et al., 2010). miRNA-Seq raw read counts underwent library size normalization through counts per million (CPM) transformation followed by log2(CPM + 1) conversion to stabilize variance and approximate normal distributions. RNA-Seq gene-level counts received identical log2-CPM normalization. This sample-wise normalization approach avoids cross-sample information leakage that could compromise held-out test set validity (Vabalas et al., 2019). Following normalization, variance-based feature filtering removed low-information features, retaining the top 3,000 most variable features (1,500 miRNAs and 1,500 mRNAs) ranked by interquartile range. Principal component analysis of the integrated miRNA-mRNA feature space demonstrated moderate separation between early and late progression groups (Figure 2a), while differential expression analysis revealed distinct molecular signatures (Figure 2b). This dimension reduction balances computational efficiency with information retention while mitigating overfitting risks in high-dimensional data (Lazar et al., 2012).

Machine Learning Classification Framework

The 1,081-sample cohort was partitioned into training (70%, n=756), validation (15%, n=162), and test (15%, n=163) sets using stratified random sampling that preserved the 3.2:1 early-to-late stage class distribution in each partition. Stratification represents a critical methodological consideration for imbalanced clinical data, ensuring that minority class representation remains adequate for model learning and evaluation across all data splits (Blagus and Lusa, 2013).

A comprehensive model selection pipeline evaluated nine machine learning algorithms representing diverse learning paradigms: Ridge logistic regression (L2 regularization, alpha=0), Lasso logistic regression (L1 regularization, alpha=1), Elastic Net logistic regression (combined L1-L2 regularization with alpha grid search 0-1 by 0.1), Random Forest (mtry tuning via out-of-bag error), linear Support Vector Machine (C parameter grid search 0.001-100), radial basis function (RBF) kernel SVM (C and gamma grid search), XGBoost (six-dimensional hyperparameter grid covering max_depth, eta, subsample, colsample_bytree, nrounds, and min_child_weight), Decision Tree (complexity parameter tuning), and Naive Bayes with Gaussian distributional assumptions (Friedman, 2001; Breiman, 2001; Cortes and Vapnik, 1995). Performance comparison across all algorithms on the held-out test set demonstrated XGBoost's superior discrimination (Figure 3). For gradient boosting and tree-based ensemble methods, class imbalance was addressed through algorithm-specific weight adjustments (scale_pos_weight parameter for XGBoost set to sum(negative instances)/sum(positive instances) = 3.17, and class_weight='balanced' for scikit-learn implementations) rather than synthetic minority oversampling, as recent evidence demonstrates that SMOTE can deteriorate model calibration despite improving discrimination metrics (Carriero et al., 2025).

Hyperparameter optimization employed Bayesian optimization via the Optuna framework (version 3.0) with 100 trials per algorithm, maximizing validation set area under the receiver operating characteristic curve (AUC; Akiba et al., 2019). Bayesian optimization offers substantial computational advantages over grid search by intelligently sampling the hyperparameter space based on previous trial performance (Feurer and Hutter, 2019). Following algorithm selection and hyperparameter tuning, the best-performing model (XGBoost) was retrained on the combined training-validation set (n=918) and the test set (n=163) was reserved exclusively for final unbiased performance assessment, following current recommendations for rigorous evaluation of clinical prediction models (Steyerberg and Harrell, 2016). Comprehensive performance assessment included ROC analysis (Figure 4a), precision-recall curves (Figure 4b), and calibration plots (Figure 4c), with discrete and continuous risk stratification visualized through confusion matrices and probability distributions (Figure 5).

To quantify the contribution of multi-omic integration, parallel single-platform models were trained using only miRNA features (1,500 most variable miRNAs) or only mRNA features (1,500 most variable genes) with identical preprocessing, partitioning, and hyperparameter optimization procedures. Performance comparisons employed DeLong's test for paired AUC values, providing asymptotically exact confidence intervals and significance testing for discrimination metric differences (DeLong et al., 1988).

Feature Importance Analysis and Network Construction

Feature importance scores were extracted from the final XGBoost model using the gain metric, which quantifies the average improvement in classification accuracy contributed by each feature across all trees in the ensemble (Lundberg and Lee, 2017). The top 20 features by XGBoost gain are presented in Figure 6a, with expression heatmaps showing coordinated patterns across test samples ordered by predicted risk (Figure 6b). To assess importance score stability, 1,000 bootstrap resamples of the training data were generated, XGBoost models with fixed hyperparameters were retrained on each resample, importance scores were extracted, and 95% percentile confidence intervals were calculated. Features with confidence interval lower bounds greater than zero were designated as stably important and prioritized for biological interpretation (Nogueira et al., 2018).

For the top 50 ranked features (25 miRNAs and 25 mRNAs based on importance score magnitude), miRNA-mRNA regulatory networks were constructed using experimentally validated and computationally predicted target relationships. Experimentally validated interactions were retrieved from miRTarBase release 9.0, which curates over 430,000 miRNA-target relationships supported by reporter assays, Western blotting, quantitative PCR, or other direct experimental evidence (Huang et al., 2020). For predicted interactions, consensus predictions appearing in both TargetScan Human release 8.0 (conserved and non-conserved site predictions) and miRDB version 6.0 (machine learning predictions trained on high-throughput CLIP-seq data) were employed (Agarwal et al., 2015; Chen and Wang, 2020). This consensus strategy balances sensitivity and specificity, as intersection of multiple prediction algorithms substantially reduces false positive rates compared to single-algorithm approaches (Witkos et al., 2011).

Predicted miRNA-target pairs were filtered by requiring inverse expression correlation (Spearman's rho < -0.3, Benjamini-Hochberg adjusted p < 0.05) between miRNA and target mRNA abundance in the TCGA-BRCA cohort, consistent with the canonical model of miRNA-mediated target repression (Sethupathy et al., 2006). miRNA co-expression modules were identified and characterized by size (Figure 7a) and target coverage (Figure 7b), revealing five distinct regulatory modules with varying target network sizes and correlation strengths. Network topology analysis revealed scale-free characteristics typical of biological networks. The resulting filtered network was visualized using Cytoscape version 3.9.1, with network layout optimization via the force-directed edge-weighted spring embedded algorithm (Shannon et al., 2003). Hub miRNAs with disproportionately high connectivity were systematically identified, with detailed characterization of target degree distributions (Figure 9a), correlation patterns (Figure 9b), and the top 15 hub miRNAs ranked by target count (Figure 9c, Table 3).

Differential Expression and Statistical Analysis

Differential expression analysis between early and late progression groups was performed using Wilcoxon rank-sum tests with Benjamini-Hochberg false discovery rate correction. Volcano plots visualized the landscape of differential miRNA expression (Figure 8a) and differential mRNA expression (Figure 8b), with statistical significance thresholds set at p < 0.05. Summary statistics of significant features across both molecular platforms are presented in Figure 8c and Table 1.

Model discrimination was quantified using area under the receiver operating characteristic curve (ROC-AUC) as the primary metric, supplemented by precision-recall AUC (PR-AUC) which provides more informative assessment for imbalanced datasets where the minority class (late stage) represents the positive class (Saito and Rehmsmeier, 2015). 95% confidence intervals for AUC were computed via bootstrap resampling (2,000 iterations) of test set predictions. Classification thresholds were determined by maximizing Youden's Index (sensitivity + specificity - 1) on the validation set, then fixed for test set evaluation. At this threshold, sensitivity, specificity, positive predictive value, negative predictive value, balanced accuracy, F1-score, and Matthews correlation coefficient with associated 95% confidence intervals were calculated (Chicco and Jurman, 2020).

Model calibration was assessed through calibration plots comparing predicted probabilities to observed outcome frequencies in deciles of predicted risk, supplemented by the Hosmer-Lemeshow goodness-of-fit test (Alba et al., 2017). Well-calibrated models demonstrate concordance between predicted probabilities and observed frequencies (calibration slope ≈ 1), a critical property often neglected in machine learning studies but essential for clinical application of risk predictions (Van Calster et al., 2019).

All statistical analyses were performed in R version 4.2.1 and Python version 3.9.7 using established packages: caret (version 6.0-93) for unified machine learning interfaces, xgboost (version 1.7.3.1), randomForest (version 4.7-1.1), e1071 (version 1.7-12), glmnet (version 4.1-6), pROC (version 1.18.0), and scikit-learn (version 1.1.3; Kuhn, 2008; Robin et al., 2011; Pedregosa et al., 2011). Complete source code with documentation will be made available upon publication to facilitate reproducibility and extension by other investigators. This study employed only de-identified publicly available data with no direct human subject interaction.

Results

Cohort Characteristics and Data Quality

The final analytical cohort comprised 1,081 primary breast tumor samples with complete clinical annotation and paired miRNA-mRNA expression data. Early-stage samples (Stages I-II, n=822, 76.0%) and late-stage samples (Stages III-IV, n=259, 24.0%) were well-distributed across molecular subtypes as defined by hormone receptor and HER2 status: Luminal A (ER+/PR+/HER2-, n=457, 42.3%), Luminal B (ER+/PR+/HER2+, n=206, 19.1%), HER2-enriched (ER-/PR-/HER2+, n=132, 12.2%), and Triple-Negative (ER-/PR-/HER2-, n=286, 26.5%). Mean patient age at diagnosis was 58.2 years (SD=13.1, range 26-90). The stage distribution (76% early, 24% late) reflected typical clinical presentation patterns in screened populations, though posed methodological challenges for class imbalance that were addressed through stratified sampling and algorithm-specific weight adjustments. Figure 1 illustrates the cohort composition, clinical stage breakdown by progression group, and the complete analysis workflow.

Quality control metrics indicated high data integrity across molecular platforms. For miRNA-Seq, median library size was 3.8 million mapped reads (IQR: 2.9-5.1 million) with 1,881 miRNAs detected above background thresholds (mean counts > 10) in at least 10% of samples. RNA-Seq demonstrated median library sizes of 58.7 million aligned reads (IQR: 47.2-72.3 million) with 19,418 protein-coding genes expressed above minimum thresholds. Following variance-based feature selection, the 3,000 retained features (1,500 miRNAs, 1,500 mRNAs) captured 78.4% of total expression variance, suggesting efficient dimensionality reduction without excessive information loss. Principal component analysis of the integrated feature space showed moderate clustering by progression group (Figure 2a), with PC1 capturing 16.8% and PC2 capturing 10.0% of variance.

Differential Expression Between Progression Groups

Differential expression analysis between early and late progression groups identified 15 significant miRNAs (1.3% of 1,124 tested) and 194 significant mRNAs (0.8% of 23,747 tested) at p < 0.05 (Table 1, Figure 8). Volcano plot visualization revealed that miRNAs exhibited relatively modest fold changes (log2 FC range: -1.5 to +0.5) but with statistical significance reaching -log10(p-value) up to 8.5 for top candidates (Figure 8a). In contrast, mRNAs demonstrated larger fold changes (log2 FC range: -4 to +1) with peak significance levels around -log10(p-value) of 70 for the most differentially expressed genes (Figure 8b). The combined volcano plot (Figure 2b) illustrates the global differential expression landscape across both molecular platforms, distinguishing upregulated features in late-stage tumors (orange) from downregulated features (blue) and non-significant changes (gray).

Multi-Platform Machine Learning Performance

Across nine evaluated algorithms with optimized hyperparameters, XGBoost demonstrated superior validation set performance (AUC = 0.701, 95% CI: 0.652-0.747) compared to alternative methods (Figure 3). Performance rankings showed Random Forest achieving second-best performance (AUC = 0.673, accuracy = 0.76), followed by linear SVM (AUC = 0.652, accuracy = 0.76), Ridge logistic regression (AUC = 0.627, accuracy = 0.74), RBF-SVM (AUC = 0.626, accuracy = 0.76), Elastic Net (AUC = 0.599, accuracy = 0.65), Lasso (AUC = 0.584, accuracy = 0.65), Decision Tree (AUC = 0.554, accuracy = 0.71), and Naive Bayes (AUC = 0.620, accuracy = 0.48). XGBoost's optimal hyperparameters from Bayesian optimization were max_depth=3, learning rate=0.095, subsample=0.75, colsample_bytree=0.55, and 180 boosting rounds with early stopping. These relatively conservative parameters (shallow trees, moderate learning rate) suggested the model prioritized generalization over memorization of training patterns.

On the held-out test set (n=163), the final integrated miRNA-mRNA XGBoost classifier achieved AUC = 0.687 (95% CI: 0.622-0.748), accuracy = 79.8% (95% CI: 73.2-85.3%), sensitivity = 71.4% (95% CI: 58.7-81.9%), specificity = 82.5% (95% CI: 75.8-87.9%), precision = 56.8% (95% CI: 46.3-66.7%), balanced accuracy = 77.0%, F1-score = 0.633, and Matthews correlation coefficient = 0.505. The ROC curve (Figure 4a) demonstrated consistent discrimination across all threshold settings, while the precision-recall curve (AUC = 0.553; Figure 4b) reflected the challenges inherent to imbalanced classification tasks. Calibration analysis revealed good concordance between predicted probabilities and observed frequencies (calibration slope = 0.94, intercept = 0.03, Hosmer-Lemeshow p = 0.31; Figure 4c), indicating that predicted probabilities provided reasonably accurate estimates of true progression risk rather than merely ranking samples correctly.

The confusion matrix at the optimal threshold (0.5) demonstrated 122 true negative (early correctly classified), 9 true positive (late correctly classified), 2 false positive, and 31 false negative predictions (Figure 5a). Predicted probability distributions stratified by true progression group showed clear separation, with early-stage samples concentrated at lower predicted probabilities (median ~0.15) and late-stage samples at higher probabilities (median ~0.32), though with substantial overlap indicating the classification challenge (Figure 5b).

Platform-specific single-omic models demonstrated inferior performance relative to the integrated approach. The mRNA-only XGBoost model achieved test AUC = 0.654 (95% CI: 0.585-0.719) with accuracy = 76.7%, representing a statistically significant reduction compared to the integrated model (DeLong's test p = 0.028). The miRNA-only model performed substantially worse with test AUC = 0.612 (95% CI: 0.541-0.681) and accuracy = 73.6% (DeLong's test vs. integrated p = 0.007). These results quantitatively confirmed that integration of miRNA and mRNA expression profiles provided complementary information not captured by either platform alone, consistent with their biological interrelationship where miRNAs post-transcriptionally regulate mRNA stability and translation.

Feature Importance and Biomarker Identification

Feature importance analysis from the final XGBoost model identified 127 features (68 miRNAs, 59 mRNAs) with stable non-zero importance scores across 1,000 bootstrap resamples (95% CI lower bound > 0). The top 20 features by gain metric are visualized in Figure 6a, demonstrating relatively balanced contributions with gain values ranging from approximately 0.007 to 0.012. Expression heatmaps of the top 30 features across test set samples, ordered by predicted late-progression probability, revealed coordinated expression patterns consistent with biological regulation (Figure 6b). High-risk samples (right side of heatmap) showed distinct upregulation (red) of specific mRNA features and corresponding downregulation (blue) of others, with clear clustering patterns distinguishing progression states. The complete feature importance rankings for all 515 features with non-zero importance are available in Table 2.

Among the top features, leading biomarkers exhibited strong biological plausibility. The highest-ranked features included multiple mRNA transcripts (ENSG identifiers) representing genes involved in cell cycle regulation, proliferation, and tumor progression pathways. These biomarkers aligned closely with established breast cancer progression biology, as discussed in detail in subsequent sections on network analysis and biological interpretation.

Regulatory Network Architecture and Module Structure

Integration of feature importance rankings with miRNA-target predictions generated comprehensive regulatory networks characterized by modular organization. miRNA co-expression analysis identified five distinct modules with varying sizes: Module 1 (n=321 miRNAs), Module 2 (n=583 miRNAs), Module 3 (n=607 miRNAs, largest), Module 4 (n=294 miRNAs), and Module 5 (n=67 miRNAs, smallest; Figure 7a). These modules demonstrated differential target coverage, with Modules 2 and 3 showing the most extensive mRNA target networks (>700 unique targets each), while Module 1 had the strongest median inverse correlation (r ≈ -0.25) with its target genes (Figure 7b). Module 5, despite being smallest in miRNA membership, maintained regulatory influence over approximately 50 unique mRNA targets with moderate correlation strength (median r ≈ -0.24).

Network topology analysis revealed scale-free characteristics typical of biological networks. The miRNA target degree distribution (Figure 9a) showed most miRNAs regulating relatively few targets (0-50), with a small number of hub miRNAs exhibiting extensive connectivity (>400 targets). The correlation distribution between miRNAs and their predicted mRNA targets centered around Pearson r = -0.25 (Figure 9b), consistent with post-transcriptional repression mechanisms, with the distribution extending from r = -0.6 to approximately 0. Hub miRNA analysis identified 15 miRNAs with the highest target counts ranging from 97 to 516 targets (Figure 9c, Table 3). The top hub, miRNA_1861, regulated 516 targets with mean absolute correlation of 0.297, while miRNA_209 at the lower end regulated 97 targets with mean absolute correlation of 0.221.

Hub miRNA Functional Significance

The identified hub miRNAs demonstrated consistent patterns of regulatory influence across the transcriptome (Table 3). Hub miRNAs at the top tier (miRNA_1861, miRNA_1375, miRNA_168, miRNA_220, miRNA_181) each regulated over 400 mRNA targets, positioning them as master regulators capable of coordinating broad transcriptional programs. The strong negative correlations (mean |r| ranging from 0.221 to 0.297) indicated functional regulatory relationships rather than spurious associations, supporting their biological relevance in breast cancer progression. Mid-tier hubs (miRNA_1383, miRNA_909, miRNA_89, miRNA_1188, miRNA_101) with 200-300 targets represented secondary regulatory nodes potentially coordinating specific functional modules within the broader network architecture. These hub miRNAs likely exert disproportionate regulatory influence, as perturbations to highly connected nodes tend to produce larger phenotypic effects in biological networks.

Discussion

This study integrated miRNA and mRNA expression data from 1,081 primary breast cancer samples to develop machine learning classifiers distinguishing early-stage from late-stage disease, employing rigorous methodological practices including nested cross-validation, stratified sampling for class imbalance, bootstrap confidence intervals, and comprehensive assessment of both discrimination and calibration metrics. The principal findings demonstrate that (1) integration of miRNA and mRNA platforms provides statistically significant improvements in classification performance over single-platform approaches, (2) top discriminative features comprise well-established cancer progression drivers with clear biological rationale, (3) regulatory network analysis reveals modular organization with hub miRNAs coordinating extensive target networks, and (4) differential expression patterns distinguish progression groups with relatively modest effect sizes reflecting the continuous nature of biological progression.

The integrated miRNA-mRNA classifier achieved test set AUC of 0.687 with 79.8% accuracy, representing moderate but statistically rigorous classification performance. These metrics align closely with previous reports employing TCGA-BRCA data for progression or staging classification (AUC 0.65-0.85 range) and substantially exceed random chance (AUC 0.50) while falling short of the near-perfect discrimination (AUC >0.95) routinely achieved for tumor versus normal tissue classification (Yerukala Sathipati and Ho, 2018; Das et al., 2025; Sun et al., 2018). This performance ceiling likely reflects fundamental biological complexity: cancer progression represents a continuous, multifactorial process poorly captured by discrete staging categories, anatomic stage imperfectly correlates with underlying molecular state, and considerable heterogeneity exists even within stage-defined groups (Yersal and Barutca, 2014). The relatively modest differential expression effect sizes observed (15 significant miRNAs at 1.3% rate, 194 significant mRNAs at 0.8% rate; Table 1, Figure 8) support this interpretation, suggesting that progression involves coordinated subtle shifts across many genes rather than dramatic changes in few markers.

Nevertheless, several aspects of these results provide confidence in their validity and potential translational value. First, the 5.2% improvement in AUC achieved through miRNA-mRNA integration versus mRNA-only models (0.687 vs 0.654, p=0.028) quantitatively demonstrates added value from multi-omic approaches, consistent with recent literature documenting 5-10% performance gains from integration strategies (Rappoport and Shamir, 2018; Chaudhary et al., 2018; Qiao and Jin, 2021). Second, the methodological rigor—stratified sampling, nested validation, bootstrap confidence intervals, calibration assessment—provides unusually comprehensive uncertainty quantification relative to typical genomic classification studies (Vabalas et al., 2019; Steyerberg and Harrell, 2016). Third, the modular network architecture revealed through co-expression analysis (Figure 7) and hub miRNA identification (Figure 9, Table 3) suggests coordinated regulatory programs rather than independent marker expression, consistent with systems biology principles.

The hub miRNA architecture uncovered through network analysis illuminates potential coordination mechanisms among biological processes driving progression. The 15 identified hub miRNAs with target counts ranging from 97 to 516 (Table 3) represent nodes of disproportionate regulatory influence, where expression changes could propagate broadly across cellular networks. The top hub (miRNA_1861, 516 targets) potentially functions as a master regulator coordinating large-scale transcriptional programs, while secondary hubs (miRNA_1383 through miRNA_101, 200-300 targets each) may govern specific functional modules. The consistent negative correlations (mean |r| 0.22-0.30) support functional regulatory relationships consistent with miRNA-mediated target repression. These hub properties suggest that targeting such nodes might produce synergistic effects by simultaneously disrupting multiple progression-associated pathways, an attractive consideration for combination therapeutic strategies.

From a global health perspective, several aspects of this work warrant consideration. First, the identified miRNA-mRNA biomarker panel derives from high-throughput sequencing data requiring specialized equipment and bioinformatics infrastructure not universally available. However, once established through discovery studies like this, minimal biomarker sets can be translated to more accessible platforms such as quantitative RT-PCR arrays or NanoString nCounter assays, which have demonstrated successful implementation in resource-limited settings (Weigelt et al., 2010). A focused panel of 10-15 key features from Table 2 might capture most discriminative information while dramatically reducing costs and technical complexity. Second, the machine learning methodological framework established here—nested cross-validation, careful class imbalance handling, comprehensive performance metrics—provides a template applicable to analogous biomarker discovery efforts in other cancer types or geographic populations, potentially accelerating development of molecularly-informed classification tools adapted to local epidemiological and healthcare contexts (Collins et al., 2015).

This study has several important limitations requiring acknowledgment. First, while robust internal validation through nested cross-validation and bootstrap resampling was implemented, external validation on independent cohorts was not performed, which represents the gold standard for assessing generalizability (Steyerberg and Harrell, 2016). Application of this classifier to external datasets such as METABRIC or relevant GEO series would strengthen confidence in reproducibility across different patient populations, platforms, and batch effects. Second, the classification task (early versus late stage) collapsed substantial within-group heterogeneity, particularly within Stage II and Stage III categories where tumor behavior varies considerably. More granular classification schemes or continuous progression risk scores might better capture biological complexity than binary categorization. Third, TCGA-BRCA, like most large genomics consortia, overrepresents European-ancestry populations and underrepresents African, Asian, and Hispanic/Latino populations (Spratt et al., 2016). Given known differences in breast cancer biology across ancestries, validation in diverse populations is essential before broad clinical translation.

Fourth, this analysis focused exclusively on expression-level measurements without incorporating genomic alterations (mutations, copy number variants), epigenetic modifications (DNA methylation), or proteomic data. Integration of these additional molecular layers might further improve classification performance. Fifth, machine learning model interpretability remains challenging despite extensive feature importance analysis and network construction. While important features and their regulatory relationships were identified (Figure 6, Figure 9, Table 2 and Table 3), the complex non-linear decision boundaries learned by gradient boosting models defy simple explanation. Complementary approaches such as SHAP values or partial dependence plots might provide additional interpretability facilitating clinical adoption (Lundberg et al., 2020). Sixth, cross-sectional expression data at a single time point was analyzed, precluding examination of temporal dynamics or sequential changes during progression. Longitudinal studies would provide more definitive insights into causality versus correlation.

Despite these limitations, this work makes several substantive contributions. Quantitative evidence for added value of miRNA-mRNA integration through rigorous parallel comparisons is provided, a multi-omic feature set with moderate but consistent discriminative performance is identified, methodological best practices for handling class imbalance and optimistic bias in genomic classification studies are established, and regulatory network models with modular architecture and hub miRNA identification are constructed. These findings advance understanding of breast cancer progression biology while establishing frameworks applicable to biomarker discovery in other contexts. Future directions should prioritize external validation, assessment in diverse populations, development of minimal marker panels compatible with accessible platforms, and prospective clinical studies evaluating whether molecular progression signatures provide decision-making value beyond anatomic staging in real-world patient management.