Four problems related to the delineation of individual breast cancers with respect to the expression data of patient samples were considered:

Is the patient sample ‘cancer’ or ‘normal’?

A generalized workflow for the problems is depicted in Figure 1.

Preprocessing was done in a manner similar to Sarathi and Palaniappan (Sarathi and Palaniappan, 2019). The source dataset for all problems modeled here was obtained from the TCGA. Normalised BRCA expression data was acquired using the firebrowse portal (Summary, 2016) (gdac.broadinstitute.org_BRCA.Merge_rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.Level_3.2016012800.0.0. tar.gz), and RSEM counts were obtained. The patient barcode was matched with the clinical data (gdac.broadinstitute.org_BRCA.Merge_Clinical.Level_1.2016012800.0.0. tar) to extract the patient. stage_event.pathologic_stage variable values that encode the AJCC TNM staging (Giuliano et al., 2018). The sub-stages were then merged to obtain the macro stage categories. Table 1 shows the distribution of sample stages for the breast cancer samples according to the AJCC staging system. It is noted that early-stage BC indicates TNM stage-I or stage-II cancer. Stage-III BC (including T3N1, T4, N2-3) represents loco-regionally advanced BC, whereas T3N0 represents a borderline diagnosis between stages II and III. For the purposes of our study, stages I, II, and III were combined into the ‘non-metastatic’ class.

The immunohistochemical (IHC) status of oestrogen receptor (ER) and progesterone receptor (PgR), human epidermal growth factor receptor 2 (HER2) oncogene, and Ki-67 (a marker of cell proliferation) are used together to subtype breast tumors into Triple-negative breast cancer (TNBC), HER2-positive, Luminal A and Luminal B (Giuliano et al., 2018; Dai et al., 2015), as shown in Table 2. Where reliable Ki-67 measurements are not available, an alternative assessment of tumor proliferation such as tumor grade could be used to distinguish between ‘Luminal A’ and ‘Luminal B’ (which tends to be HER2 negative). Complete ER, PgR and HER2 IHC metadata were available for 719 samples of the TCGA Breast Cancer dataset, and of these, no sample had information on the Ki-67 labeling index nor on the tumor grade, precluding precise differentiation of luminal subtypes of breast cancers into ‘Luminal A’ or ‘Luminal B’. The luminal subtypes A and B were perforce lumped into one ‘Luminal’ type. The 719 samples were accordingly annotated as 567 ‘Luminal’ (generally Luminal A with Grade 1 or 2 and Luminal B with G3), 115 TNBC (generally Grade 3), and 37 HER2 (generally Grade 3) based on the status of ER, PgR and HER2 extracted from the clinical file (Table 2).

The two most common histological subtypes of breast cancer are infiltrating ductal carcinoma (IDC - no special type) and infiltrating lobular carcinoma (ILC) (Weigelt et al., 2010). ILC tends to be difficult to diagnose, with MR imaging required for determining size and multifocality including contralateral breast (mirror image), and preferential spread to gastrointestinal tract and peritoneum (Winchester et al., 1998). The sample histological subtype is encoded in the clinical metadata ‘patient.histological_type’ with the major values being, ‘infiltrating ductal carcinoma (IDC)’ and ‘infiltrating lobular carcinoma (ILC)’, and minor values including ‘mixed histology’, ‘metaplastic carcinoma’, ‘mucinous carcinoma’, ‘medullary carcinoma’, and ‘other (specify)’.

Genes that had minimal variation in expression across the samples (i.e., σ < 1) were removed. Cancer samples which were missing stage annotation details were removed. The expression dataset was subjected to variance-stabilization using voom function in limma (Law et al., 2014). Linear modeling was then performed. The resulting dataset was split 80:20 into a training set and a holdout testset stratified on the outcome variable of each problem. It is noted that the training dataset for Problem #2 suffered an imbalance in the distribution of the outcome classes (16 metastatic vs. 837 non-metastatic samples), which prompted the application of SMOTE correction (Chawla et al., 2002) (Synthetic Minority Oversampling TEchnique; with arguments: perc. over-represented = 1,000% and perc. under-represented = 300%). Data preprocessing and analysis was done using R (www.r-project.org). The annotated pre-processed final dataset is available as Supplementary File S1.

Feature spaces for each problem were constructed using only the training dataset. Initially the differential expression of genes across cancer stages relative to healthy samples was studied using linear modelling with limma (Ritchie et al., 2015): y = α + β 1 X 1 + β 2 X 2 + β 3 X 3 + β 4 X 4 (1)

Where the independent variables are indicator variables of the sample’s stage, the intercept α is the baseline expression estimated from the controls, and β i are the estimated stagewise log fold-change (lfc) coefficients relative to controls.

We then applied a two-level contrast protocol (Muthamilselvan et al., 2023), viz. level-I: stage vs. control and level-II: inter-stages contrast, to produce the following classes of features:

Stage-salient genes obtained from all possible pairwise contrasts between the cancer stages using the following model:

Where the controls themselves constitute one of the indicator variables ( X 0 ), and the β i are coefficients estimated from samples of the corresponding annotation only.

Monotonically expressed genes obtained from strictly increasing or strictly decreasing mean expression across the cancer stages.

In addition, expression contrasts specific to the problem under consideration were used, namely:

contrast of non-metastatic vs. metastatic cancers using the following model modified from Equation 2:

Where the μ i are coefficients estimated from samples of the corresponding annotation only.

three-way pairwise contrasts between the molecular subtypes; viz. (i) Luminal vs. HER2+, (ii) Luminal vs. TNBC and (iii) HER2+ vs. TNBC using the following model modified from Equation 2:

Where the δ i are coefficients estimated from samples of the corresponding annotation only.

contrast of ductal vs. lobular histologies using the following model modified from Equation 2:

Where the ϑ i are coefficients estimated from samples of the corresponding annotation only.

The above strategies yielded problem-specific chimeric feature spaces that could span the informative dimensions in each case.

A composite feature space comprising the top-ranked genes from the linear model, stage-salient genes, and genes from the problem-specific contrast was subjected to the consensus of two feature selection techniques: (i) Boruta, a wrapper algorithm using Random Forest to select features based on a measure of importance to the outcome variable of interest (Kursa and Rudnicki, 2010); and (ii) Recursive Feature Elimination (RFE), a method that uses backward selection passes to trim the space of predictor variables. The workflow of the machine learning model development in Figure 1 presented in the context of cancer v/s normal was adapted for the non-metastatic v/s metastatic, molecular subtype, and histological subtype classification problems. The training dataset with the final set of features was loaded onto models based on six different algorithms, including Random Forest (ensemble bagging classifier that builds numerous decision trees and ‘bags’ the majority vote), Support Vector Machine (geometric method that finds the maximum margin separating hyperplane in high-dimensional space), k-NN (based on distance-based proximal classes), 1-layer and 2-layer Neural Networks, and XGBoost (ensemble boosting classifier that builds a sequence of classifiers iteratively ‘boosted’ on challenging instances).

Subsequent to an 80:20 train-test split, algorithm-specific hyperparameter configuration was optimized using 10-fold cross-validation on the training dataset for each of the six algorithms considered. Different algorithm classes were then compared based on their outer-fold testset performance, to identify the optimal algorithm class for each learning problem. The design of such a nested model selection prevents information leakage between model tuning and evaluation, and provides for a more reliable assessment of model generalizability to unseen cohorts than merely cross-validation. Evaluation metrics on the holdout testset as well as external datasets (described below) included balanced accuracy, F1-score, area under ROC (AUROC), Mathews’ correlation coefficient (MCC), and Positive Predictive Value (PPV).

The overall best model for each problem was validated primarily by performing inference on out-of-domain external datasets. Table 3 shows the datasets used in the development and validation of the ML models for the respective classification problems. In addition, we sought to obtain concurrence for our models from multi-omic signatures, as discussed below.

A prediction pipeline that integrates the predictions from all the models into one combined readout was designed. A schematic for one such cascade model is shown in Figure 2. Based on the decision at the shown fork, the new sample may be taken forward for assessment of metastatic potential and molecular/ histological subtyping. The final readout for a sample from the cascade classifier would consolidate the inference from each model; for e.g., ‘Metastatic triple-negative ductal cancer’. This formed the basis for the development of BC-Predict.

The workflow for this learning problem is shown in Figure 1. Stratified sampling of the TCGA BRCA dataset based on the class ‘cancer’ or ‘normal’ yielded a training dataset of 90 Normal and 854 Cancer samples, and a test dataset of 22 Normal and 212 Cancer samples. The 24 stage-salient genes from the contrasts shown in Equation 2 (namely CHRNA6, MMP10, DEPDC1, COX7A1, KCNK15, MFSD4, CDH19, CXCL5, AKR7A3, DEGS2, CST2, LOC100124692, GDF5, FOXA1, EGR3, FOS, FOSB, DUSP1, FREM1, EGR1, HFM1, ABCA10, KLK5, KCNA1) were combined with the top 10 linear modelling genes from Equation 1 (namely NEK2, MMP11, PKMYT1, GPAM, CPA1, COL10A1, MYOC, KIF4A, CA4, LYVE1) to obtain 34 base features for feature selection. Application of the RFE procedure identified ten features for model development, including two stage-salient genes (FREM1, ABCA10) and eight genes from the linear model (NEK2, MMP11, PKMYT1, GPAM, CPA1, COL10A1, CA4, LYVE1). Of the six ML models trained, four models yielded >99% balanced accuracy on the training set. Subsequent evaluation on holdout testset identified only one model class with 100% accuracy, namely the neural network with one hidden layer model (Supplementary File S8). The model was re-built using the full dataset and validated on external datasets: (i) BRCA-KR, yielding a balanced accuracy ∼94.00%; and (ii) GTEx, yielding ∼100% accuracy (all correct predictions). Together, the model yielded an overall balanced accuracy ∼97.42% on external validation (Table 7). The details could be found in Supplementary File S9, along with the prediction probabilities for all instances in both the external validation. Prediction probability is a measure of the strength of evidence for the predicted class, and based on the distribution of its values, recommendations for evidence of the predicted class may be generated. It was observed that correct predictions were supported by very strong prediction probabilities (>0.9) relative to incorrect predictions.

The workflow for this learning problem is a variation on Figure 1, and available in Supplementary File S10. Stratified sampling of the TCGA BRCA dataset based on the class ‘non-metastatic’ or ‘metastatic’ yielded a training dataset of 837 non-metastatic and 16 Metastatic samples, and a test dataset of 209 non-metastatic and 4 Metastatic samples. SMOTE balancing of the training dataset yielded a dataset with 480 non-metastatic and 176 Metastatic samples. The contrast shown in Equation 3 between non-metastatic and metastatic samples in the SMOTE-balanced dataset produced two lists of genes, one sorted by log-fold change and the other by significance (adj. p-value). The consensus of the top 50 genes from the two lists identified 15 features (namely SRMS, OXT, MMP27, LOC158696, C4orf26, CECR4, ANKRD55, GALNTL6, KRTAP3-1, FAM69C, AFP, CCDC33, SLC5A5, CXorf48, RGS7), to which were added the six top genes by significance missing in the consensus (namely GIP, SSX5, LOC100101938, C9, ASZ1, COX8C). Finally, these 21 genes were pooled with the 24 Stage-salient genes discussed in Cancer V/s Normal classification problem, to obtain 45 base features for feature selection. Application of the Boruta protocol yielded 14 features, while application of RFE procedure yielded just five features. The five RFE features were a subset of the features identified by Boruta, thus we obtained five consensus features for model development, namely DEPDC1, FOSB, DUSP1, MMP27 and ABCA10. Of the six different ML models trained, three models yielded >99% balanced accuracy on the training set. Subsequent evaluation on the holdout testset identified the neural network with one hidden layer model as the best performing model class, with 82.24% balanced accuracy (Supplementary File S8). The model was re-built using the full dataset and validated on the BRCA-KR and GSE18549 datasets, yielding an overall balanced accuracy ∼88.22% on the external validation (Table 7). The details could be found in Supplementary File S10, which includes the prediction probabilities for all instances in the external validation. On inspection of the distribution of prediction probabilities, correct predictions were found to be supported by high values (>0.75) relative to incorrect predictions.

The workflow for this learning problem is a variation on Figure 1, and available in Supplementary File S11. Stratified sampling of the TCGA BRCA dataset based on the molecular subtype class (‘Luminal’ or ‘TNBC’ or ‘HER2’) yielded a training dataset of 434 Luminal, 30 HER2 and 92 TNBC samples, and a test dataset of 113 Luminal, 7 HER2 and 23 TNBC samples. The three-way pairwise contrasts shown in Equation 4 between the molecular subtypes; viz. (i) Luminal vs. HER2, (ii) Luminal vs. TNBC and (iii) HER2 vs. TNBC; yielded subtype-specific genes, from which the top ten genes of each subtype (by significance) were pooled together to obtain 30 base features for feature selection (namely MLPH, AGR3, CA12, TBC1D9, AGR2, TFF3, SIDT1, FZD9, BCAS1, CXorf61, ERBB2, PGAP3, STARD3, C17orf37, GRB7, PSMD3, PCSK6, PNMT, TCAP, LOC150622, GATA3, ANXA9, FLJ45983, PRR15, FOXA1, DEGS2, SLC44A4, ZMYND10, KCNK15, NAT1). Application of the Boruta protocol did not identify any redundant feature, whereas application of RFE procedure yielded 16 features. These 16 features were identified as the consensus features for model development, namely GATA3, AGR3, CA12, TBC1D9, ERBB2, MLPH, KCNK15, ANXA9, FLJ45983, GRB7, PGAP3, STARD3, SLC44A4, PCSK6, FOXA1 and BCAS1. Of the six different ML models trained, the Random forest model provided superlative performance on both the training and outerfold test sets, with balanced accuracies of >99% and 91.43% respectively (Supplementary File S8). The model was re-built using the full dataset and was validated on the METABRIC dataset, yielding a balanced accuracy ∼88.79% (Table 7). Availability of the TNBC-only dataset provided an opportunity to execute a second out-of-cohort validation, yielding correct identification of 25 TNBC samples out of the total 26 samples (96.15% accuracy). The details could be found in the Supplementary File 11, including the prediction probabilities for all instances in the METABRIC and TNBC external validation datasets. On inspection of the distribution of prediction probabilities, correct predictions were found to be supported by high values (>0.7) relative to incorrect predictions. We investigated the 16 features used in the RandomForest model for feature importance based on mean decrease in Gini score in R caret (Kuhn, 2008). The top five features contributing to the model performance were identified as GATA3, CA12, AGR3, TBC1D9, and MLPH (Figure 5).

Stratified sampling of the TCGA BRCA dataset based on the histological subtype (‘IDC’ or ‘ILC’) yielded a training dataset of 624 IDC and 162 ILC samples, and a test dataset of 156 IDC and 40 ILC samples. The contrast shown in Equation 5 between the ductal and lobular histologies was used to detect differentially expressed genes between the two histologies, specifically applying a log-fold change threshold, |lfc| >2, to binarize genes useful as features. This obtained 62 base features for feature selection. Application of the Boruta protocol yielded 58 features, while application of the RFE procedure yielded 24 features. The 24 RFE features were a subset of the Boruta features, thus we obtained 24 consensus features features for model development, namely ADCY5, ALDH1L1, ANKRD43, C1orf64, C7, CAPN8, CCL14, CDH1, CIDEA, CTSG, DARC, F7, FXYD1, HPX, IGFN1, MMP1, PEBP4, PLCXD3, PROL1, SHROOM1, TFAP2B, TFF1, TNNT3, and WNK4. Of the six different ML models trained, four models yielded >95% balanced accuracy on the training set. Subsequent evaluation on the holdout testset identified XGBoost as the best performing model class, with 84.94% balanced accuracy (Supplementary File S8). To mitigate overfitting to the larger IDC class at the expense of the ILC class, we sought to combine the XGBoost model with the 1-layer neural network model, producing a voting ensemble classifier with a slightly better 88.74% balanced accuracy on the holdout testset (Table 7). The ensemble model re-built using the full dataset was validated on the external dataset: brca_mbcproject_wagle_2017, encoding both the histological subtypes of interest (IDC and ILC) as well as other subtypes such as ‘mixed histology’, ‘DCIS’ (ductal carcinoma in situ), and ‘NOS’. Predictions were accepted if the two models of the ensemble agreed on the predicted class. If the models disagreed on the predicted class, then the predictions were rejected as ambiguous. Such instances represent challenges to the ensemble classifier whose resolution might not be simple. Omitting the eleven such instances from the external dataset, we obtained correct predictions on all 91 IDC samples as well seven (out of thirteen) ILC samples, yielding an ensemble accuracy ∼94.23% and balanced accuracy ∼76.92% (Table 7). Even with ensembling, generalization errors persisted in learning the ILC class, with an imbalance in the type-II error between the two classes. The details could be found in Supplementary File 12, including the prediction probabilities for all instances in the external validation. On inspection of the distribution of prediction probabilities, correct predictions were found to be supported by high values (>0.7) relative to incorrect predictions. Histological subtyping from molecular features has remained a refractory learning problem, and we have made our models and code freely available for non-commercial use (www.github.com/apalania/BC-Predict_Histological).

Stage-salient miRNA were identified using the two-level contrasts of the miRNA expression data, and then their targets were identified using the R multiMiR library (Supplementary File S13). Based on these results, we determined the concordance between the regulatory miRNAs and their target genes. Temporal concordance in expression exists if the salience in miRNA expression is at least as early as the salience in target gene expression. If the expression pattern of miRNA is discordant with its target gene, a paradoxical aberration with a protective function is possible. Table 8 summarizes the validation of stage-salient gene expression from the angle of miRNA expression. Concordance between the mRNA and miRNA in the direction of expression as well as the temporal dimension is achieved for 13 stage-salient genes: MMP10, DEPDC1, CDH19, FOXA1, DEGS2, CST2, AKR7A3, EGR1, EGR3, FOS, FOSB, FGF2, and HCN2. The key regulatory miRNAs decoded by stage included 25 stage-salient miRNAs (Supplementary File S13), appearing to regulate most of the stage-salient genes. Stage-salient miRNA that were fully concordant with target mRNAs included hsa-miR-182-5p, hsa-miR-210-3p, hsa-miR10b-5p, hsa-miR-200a-5p, hsa-miR-96-5p, hsa-miR-21-5p, hsa-miR-133a-3p, hsa-miR-335-5p, hsa-miR-204-5p, and hsa-miR-145-5p. Further, four of the stage-salient miRNAs regulated genes that featured in the ML models, namely hsa-miR-210-3p, hsa-miR10b-5p, hsa-miR-200a-5p, and hsa-miR-96-5p. Only five stage-salient miRNAs displayed no overlap between their targets and stage-salient genes, and conversely, eleven stage-salient genes were predicted to be free of regulation by a stage-salient miRNA (namely COX7A1, DACT2, KCNK15, MFSD4, DSC3, KLK5, KRT15, LOC100124692, ABCA10, MAPK8IP2, and MASP1). The complete and fully detailed analysis could be found in Supplementary File S13.

Aberrant methylation in the core/ proximal promoter regions as well as enhancers could have profound regulatory effects on gene expression. We obtained a total of 22 stage-salient DMGs from the consensus of Averep, Mvalue, and MethylMix procedures: 1 stage-I salient DMG (VOPP1), 8 stage-II salient DMGs (HS3ST3B1, CPLX1, EGR1, GMDS, ITPKB, TGFB1I1, C6orf145, SHC1), 10 stage-III salient DMGs (BTLA, TNFAIP2, PHYHIPL, LYN, MAML2, C16orf62, GPRC5B, CAPN9, AIPL1, AGAP1), and 4 stage-IV salient DMGs (CNP, TSPYL5, SLC7A5, HCN2). Salient methylation of a gene is an epigenetic mechanism to tune gene expression and would precede changes in its expression. In this respect, the stage-II salient methylation of EGR1 possibly set the stage for its stage-IV salience (minimization) in expression. It is observed that the stage-IV salient hypermethylation of HCN2 was at odds with its stage-IV salient overexpression.

Mining the methylation patterns of all stage-salient genes for differential methylation-driven genes revealed five transcriptionally predictive genes negatively correlated with gene expression, namely AKR7A3, COX7A1, DEGS2, EGR1, and FOXA1 (Figure 6). Four of these genes exhibited two-component mixtures of methylation distribution, indicating a probable shift in methylation levels in cancer samples relative to healthy ones. COX7A1 showed three-component mixtures of methylation distribution, indicating a reliance on methylation to achieve regulatory fine-tuning. Table 9 summarizes the methylation patterns for these five genes, showing the correlation size with expression and if the correlation is concordant as well. In the epigenetic context, the methylation pattern of a gene could be deemed concordant with its expression if maximal methylation is observed ahead of minimal mRNA expression. FOXA1 mRNA expression is at odds with both its epigenetic profiles (methylation and miRNA), suggesting that epigenetic modulation was being used to restore FOXA1 aberrant expression. Concordance in methylation is observed for AKR7A3, DEGS2, EGR1, and COX7A1, providing strong support for their stage-salience. The above genes except COX7A1 were also concordantly modulated by stage-salient miRNAs. Such findings lead to a belief in the existence of concert between the different layers of omics, adding ‘definiteness’ to gene expression on the path to phenotypic states. Further investigations could shed light on the emergent hypotheses in the future. The mixture decomposition of methylation patterns of the remaining stage-salient genes is provided in Supplementary File S14. It could be seen, for e.g., that the methylation of ABCA10 is positively correlated with its expression, escaping clear interpretation.

We searched Pubmed (www.pubmed.gov) using the keyword: “breast cancer” AND “stage specific” AND “gene”, and found a handful of known stage-specific genes. TIEG (or KLF10) is an anti-metastasis/ tumor-suppressor gene, which inhibits invasive breast cancer by blocking EGFR transcription in the EGFR signalling pathway (W et al., 2012). Stage-specific expression of KLF10 in breast cancer biopsies has been published, with sustained downregulation leading to complete absence of expression in invasive subtypes (Subramaniam et al., 1998). Here KLF10 expression is found to be decreasing with stage relative to the normals. γ-Synuclein (SNCG) expression is strongly correlated with the stages of breast cancer, showing little expression in normal or benign samples and increasing expression with cancer stage, and detectable only in a subset of patients (Wu et al., 2003). Here we find increasing expression of SNCG in late-stage cancers, but downregulated with respect to expression in normal samples, which is a contrarian finding.

The distinction between IDC and ILC has previously frustrated learning algorithms. An XGBoost model with 147 clinical, histopathological, mammogram features, and sonographic features has been reported with an internal testset accuracy of 0.84 on the binary classification problem (Vy et al., 2022). An AutoML deep-learning approach for identifying IDC samples alone from whole slide images yielded 0.85 accuracy on an independent dataset (Zeng and Zhang, 2020). Another study for classifying IDCs as early-stage vs. late-stage yielded an AUROC of 0.47 on the external validation (Roy et al., 2020). In this context, the external validation of our model yields a significant improvement on the state-of-the-art. However the limited sensitivity to ILC samples (conversely, specificity to IDC samples) in the external dataset presents an outstanding challenge in the histological classification of breast cancer from molecular information. Some noteworthy features from this model include: (i) CDH1 (E-cadherin), whose germline mutations were strongly associated with lobular carcinoma (Corso et al., 2018), was found to have a specific downregulated expression signature in ILC samples; (ii) CCL14, which is known to promote angiogenesis and metastasis in breast cancer (Li et al., 2011), was found oncogenic in expression across both histological subtypes. Further improvements to histological subtyping models could come from:

stacking the classifiers: the ensemble of XGBoost and neural network used herein showed that the classifiers disagree on many instances preventing a consensus classification. In such cases, improvements to the performance tradeoff could be achieved by ‘weighting’ the contribution of the two constituent models to the final prediction.

Available genomic assays (commercial or otherwise) for prognosticating breast-cancer adjuvant chemotherapy include the following gene-signature panels:

Prosigna (50 genes from PAM50 for intrinsic subtype classification, 8 housekeeping genes used for signal normalisation, 6 positive controls, and 8 negative controls)

Scanning the signatures in these genomic assays against the ten features used in our ‘normal’ vs. ‘cancer’ model yielded: two genes in common with Prosigna (FOXA1, MMP11), one gene with OncotypeDX (MMP11), one gene with HER2DX (NEK2), and one gene with Breast Cancer Index (NEK2). Scanning these signatures against the 16 features used in our molecular subtyping model yielded: four genes with Prosigna (ERBB2, FOXA1, GRB7, MLPH), four genes with HER2DX (ERBB2, GRB7, STARD3, AGR3), two with OncotypeDX (GRB7 and ERBB2), and two with Guardant360 (ERBB2, GATA3). Scanning these signatures against the 24 features used for histological subtyping yielded: one gene with Guardant360 (CDH1). Scanning these signatures against the five features used in the non-metastatic vs. metastatic model did not identify anything in common. These results indicate that the models developed in this work are novel and deserving of clinical validation. A summary of the existing gene-signature diagnostic tests (with their indications and outcomes) together with a comprehensive comparative study is provided in Supplementary file S19.

To transition the results obtained from our studies, we developed BC-Predict which serves the models developed in a cascade inference engine and provides a comprehensive characterization of the given sample (Figure 2). The BC-predict web-server is built on Rshiny (Beeley, 2016) and deployed for academic research at https://apalania.shinyapps.io/BC-Predict. All predictions are accompanied by prediction probabilities to provide confidence for the predicted class. Documentation and video tutorial for the use of BC-Predict are also provided. BC-Predict generates a unified readout that could nominally support medical decision-making contingent to clinical validation and further refinement. An alternative modeling process that used a nested stratification structure instead of sequential stratification was also investigated but did not yield an improvement. Though the cancer vs. normal model improves on the benchmark, iterative refinement and better datasets could yield further performance improvements for all models. Below we present a systematic enumeration of the limitations of our models and suggested coping strategies:

The metastatic model does not distinguish among the stages in pre-metastatic cancer. A refinement may be necessary to discriminate between the early-stage cancers (stages I and II) and stage-III cancers among the pre-metastatic cancers.

Gene-signature methods remain the clinical standard for both their effectiveness and utility, and works such as ours are a step forward in resolving difficult challenges. Such diagnostic models need to be clinically validated and approved for use by national regulatory bodies such as the FDA (Food and Drug Administration, USA), MHRA (Medicines and Healthcare products Regulatory Agency, UK), EU MDR (European Union Medical Device Regulation), NMPA (National Medical Products Administration, China) and CDSCO (Central Drugs Standard Control Organization, India). Models are complicated by cohort selection bias; for e. g., breast cancer in Black population presents in younger patients and more difficult to treat forms (aggressive, grade-III, TNBC or HER2+) than in Hispanic population, with poorer prognosis. Also, metastatic breast cancer is rarely synchronous (more metachronous) in developed nations as opposed to metastatic cancer on presentation in emerging nations. In addition to these variations, AI-based diagnostic modalities need to contend for the interplay of risk factors that could enable or confound the predictions: pre-menopausal vs. post-menopausal, node-positive or not, complete hormonal profile and NPI score. Clinical validation of BC-Predict would involve the synthesis and use of specific forward and reverse primers for each model feature to perform qRT-PCR on the isolated RNA of resected biopsy sample from a patient. Post-quantification (normalized counts) and log₂ transformation, the inference model may be served to yield a prediction. Prior to such deployment, calibration of qRT-PCR may be necessary and could involve reference genes as used in, say, NOVAprep-miR-Cervix (Kniazeva et al., 2023).

In summary, we have developed performant de novo models to characterise breast cancer heterogeneity agnostic of hypothesis. The candidate stage-salient biomarkers could play a role in the progression of breast cancer, whose varying manifestations underlie differential response to treatment regimens. Developing models from minimal feature spaces has several advantages, chief among them being sensitivity to heterogeneous individual presentation, and generalization to out-of-domain population. One example of this in the present study is the performant external validation of the Molecular Subtype model on the TNBC-only African-enriched multiethnic international cohort (25/26 samples correctly identified). It is noteworthy that TNBC is also the most common molecular subtype in the Indian subcontinent, and has frustrated drug discovery programs with few druggable targets. It may be noted that the use of mere five features in the metastatic model mitigates against the limited datasets available, and offers realistic prospects for useful generalization in clinical diagnostics. Validation analysis with miRNA strongly supported DEPDC1, FOSB and DUSP1 as potential biomarkers for metastasis. More generally, the candidate model features identified here could provide novel hypotheses for chemotherapy and immunotherapy investigations. We would like to acknowledge that the late-integration of multi-omics has not consistently provided conclusive evidence for the features used in the models, yielding possible directions for future investigations. Our study overcomes certain limitations of earlier models, namely reporting of balanced performance metrics, availability for academic research, and inclusion of external validation. The confidence returned by BC-Predict predictions could be used to safeguard against weak and uncertain evidence, addressing the hazard with AI/ML modelling (Yao et al., 2022). The clinical translation of AI/ML models would be a step forward for personalized medicine, necessitating adequate regulation to ensure the benefits of AI for all (El Naqa et al., 2023; Hickman et al., 2021). Validation and assurance of model quality could alleviate the risks of distribution drift and cohort selection bias, and pave the way for clinically effective decision support aids in precision oncology centers. The realisation of software-as-medical-devices promises to revolutionize the diagnosis, triage, and treatment of cancers.