{"paper_id":"5333d960-2a42-430c-bea4-cd142a485675","body_text":"RNA foundation models enable generalizable endometriosis disease 1 \nclassification and stable gene-level interpretation 2 \n 3 \nNiccolo McConnell1*, Jennifer Kelly1*, Ravikanth Tadikonda2, Natalia Mulligan1, Joao Bettencourt-Silva1, 4 \nMatthew Madgwick1, Ritesh Krishna1, James Strudwick1, Ashley Evans1, Stephen Checkley2, Anna Paola 5 \nCarrieri 1, Michalis Smyrnakis2, Charles H Knowles3, Laura-Jayne Gardiner1Y 6 \n 7 \n* Authors contributed equally  8 \n 9 \n1 IBM Research Europe, The Hartree Centre, Sci-Tech Daresbury, WA4 4AD, UK,  10 \n2 Hartree Centre, STFC Daresbury Lab, Daresbury WA4 4AD, UK. 11 \n3Cleveland Clinic London, 33 Grosvener Pl, London SW1X 7HY 12 \n 13 \nY Corresponding author: Dr. Laura-Jayne Gardiner (laura-jayne.gardiner@ibm.com) 14 \n 15 \nAbstract  16 \n 17 \nEndometriosis is a chronic inflammatory condition with significant diagnostic delays impacting one in ten 18 \nreproductive age women worldwide. While machine learning (ML) models trained on transcriptomic data 19 \nshow promise for disease prediction, limited generalizability across independent patient cohorts has 20 \nhindered clinical translation. Foundations models (FMs) pretrained on large-scale transcriptomic data offer 21 \npromise to learn transferrable, biologically meaningful representations that could support cross -cohort 22 \npredictions. We assembled a 12 -cohort bulk RNA -seq benchmark (334 samples) and developed a 23 \ncomputationally efficient pipeline to test whether FMs improve endometriosis classification, an approach 24 \nnot previously applied to this disease. Using AutoXAI4Omics with cohort-aware validation, we compared 25 \nembeddings derived from five state-of-the-art RNA FMs against TPM baselines. In cross-cohort prediction, 26 \nFM embeddings significantly improved performance, achieving a weighted F1 -score of 0.83 vs. 0.68 for 27 \nthe baseline. To allow gene-level interpretation of FM embedding models, we introduce classified-aligned 28 \nintegrated gradients (CA -IG), an interpretability approach aligning gene -level attributions to the 29 \ndownstream classifier without end-to-end finetuning. CA-IG revealed a conserved set of predictive genes 30 \nfrom FM embeddings across cohort-validation regimes, contrasting with unstable baseline explainability, 31 \nsuggesting that FM embeddings prioritized transferable disease-related signal over cohort-specific effects. 32 \nThese genes include novel candidates that converge on biologically plausible pathways for endometriosis. 33 \n 34 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nKeywords 35 \nFoundation Models, Machine learning, RNA, Transcriptomics, Endometriosis, disease classification, 36 \nexplainability 37 \n 38 \nBackground 39 \n 40 \nEndometriosis remains a diagnostic and therapeutic challenge, affecting one in ten women of reproductive 41 \nage and imposing a significant personal and economic burden. Definitive d iagnosis still hinges on 42 \nlaparoscopy and biopsy with expected diagnosis times of approximately 9 years in the UK  (1–3). The 43 \ndisease is characterized by ectopic endometrial -like tissue and typically presents with chronic pelvic or 44 \nabdominal pain , dysmenorrhoea, dyspareunia, fatigue, and gastrointestinal or urinary symptoms  (4).  45 \nEndometriosis is also a leading cause of subfertility/infertility, with a prevalence of 30-50% among women 46 \nseeking infertility care (5). Health-related quality of life is markedly reduced, with increased anxiety and 47 \ndepression; indirect and direct costs are considerable, and in the UK alone the annual economic burden has 48 \nbeen estimated at £8.2 billion (6). Clinical assessment is further complicated by phenotypic heterogeneity 49 \nand a poor correlation between lesion burden and symptom severity, while imaging is often inconclusive 50 \nfor superficial disease and reliable blood or tissue biomarkers remain lacking (7–9).  Hindering diagnostic 51 \ndevelopment is the fact that precise molecular drivers of endometriosis physiology are not well defined or 52 \nunderstood. As such, there is an increasing focus and need to interrogate high -throughput omics data like 53 \ntranscriptomics, which offer the ability to study genome-wide expression patterns and reveal dysregulated 54 \npathways and genes across different patient groups, lesion types and tissue/fluid types. Insights from these 55 \ntypes of studies offer potential for improved diagnostics, biomarker identification, disease subtype 56 \nclassification, increased mechanistic understanding of pathogenesis and personalized medicine approaches. 57 \nThis has motivated the application of AI to predict or model endometriosis from high-throughput expression 58 \ndata, including microarrays and RNA -seq, for various applications such as disease classification and 59 \nbiomarker discovery (9–13). 60 \n 61 \nAkter et al. used an ensemble-based approach for endometriosis prediction on single-centre bulk RNA-seq 62 \ndata (n=38) and DNA-methylation data, and reported near-perfect internal performance (F1 ~ 0.97) using 63 \nleave-one-out cross- validation (CV) (12). Revisiting the same cohort, another study used elastic-net feature 64 \nselection with boosted and bagged trees and achieved performance around F1~0.86 (5 -fold CV) (11). On 65 \nmicroarrays, Zhang et al. (2023) proposed a five-gene panel (FOS, EPHX1, DLGAP5, PCSK5, ADAT1) for 66 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nendometriosis classification and achieved a AUC of 0.836 on a held-out test set and > 0.78 on two external 67 \nvalidations (10).  A parallel line of work has pursued circulating microRNAs as non-invasive liquid-biopsy 68 \nmarkers of endometriosis, showing a six -miRNA panel exhibited a high accuracy for endometriosis 69 \nprediction (AUC=0.94) in an independent 48 -patient validation cohort (9). More recently, Bendifallah et 70 \nal. profiled the circulating miRNome in a prospective 200 -patient study and derived an 86 -miRNA 71 \nsignature with excellent internal discrimination (AUC~0.98 under cross -validation), but without multi -72 \ncentre external replication (13). As the authors note, the broader blood -biomarker literature has been 73 \nconstrained by small cohorts and methodological heterogeneity, yielding few stable, clinically adopted 74 \nmarkers. This finding was enforced by Shi et al. who demonstrated that subtype-specific tissue markers for 75 \ndeep infiltrating endometriosis were sensitive to cohort composition (14). 76 \n 77 \nWhile these studies revealed critical biological insights and significantly advanced endometriosis disease 78 \nunderstanding, their reliance on single cohorts has limited the generalizability of their models and 79 \nassociated findings to unseen datasets or cohorts. Moreover, this issue limits the clinical applicability and 80 \nvalue of their findings in real -world diagnostic settings. As such, there is a need to develop methods or 81 \napproaches which can derive more generalized insights for endometriosis that will be applicable across a 82 \nwider range of cohorts and can underpin development of robust models and derivative disease biomarkers. 83 \nA potential solution to this limitation may be the use of foundation models (FMs), large neural encoders 84 \npretrained on vast unlabeled corpora via self -supervision to learn general -purpose representations that 85 \ntransfer across tasks. Recent work has demonstrated the potential of FMs to learn representations that are 86 \nhighly transferrable across different domains, datasets and tasks. In transcriptomics, RNA FMs such as 87 \nGeneformer (15), scGPT (16),  scFoundation (17), and BulkRNABERT (18) are trained on millions of bulk 88 \nor single-cell gene expression profiles and have shown utility for downstream problems including cell-type 89 \nannotation, perturbation prediction, and disease classification . Relative to models trained from scratch, 90 \nthese encoders typically yield increased downstream task performance and greater robustness under 91 \ndistribution shift by providing richer, structure -aware features. This study represents the first systematic 92 \nevaluation of RNA foundation models for endometriosis prediction. 93 \n 94 \nIn this work we aimed to investigate if use of RNA FMs meaningfully improved generalizability of the 95 \nresulting models, allowing accurate cross-cohort disease prediction. To do this we developed a sophisticated 96 \ncross-dataset evaluation pipeline to test the domain generalizability of FM embeddings using five state-of-97 \nthe-art RNA FMs; Geneformer (15), scFoundation (17), scGPT (16), BulkRNABERT (18) and BMFM-98 \nRNA (19) and a frozen-encoder paradigm (pretrained RNA encoders are kept fixed and used to compute 99 \nembeddings which were then used as input features for a classical ML task).  The resultant pipeline is 100 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \ndetailed in Figure 1 and consisted of 1) identification of Gene Expression Omnibus (GEO) ids for relevant 101 \nbulk RNA -seq endometriosis/control datasets and automated download, extraction, amalgamation and 102 \nnormalization of transcriptomic datasets given these pre-defined ids, and LLM-powered metadata curation 103 \nwith human validation to ensure consistent and high -quality target label generation; 2) extraction of 104 \nembeddings from pre-trained frozen RNA FMs via one forward pass through each model; 3) generation of 105 \nrobust cohort-aware multi-fold strategy to evaluate within vs. cross cohort learning capabilities of using 106 \nFM embeddings as features compared to baseline normalized (TPM) gene expression counts; 4) an 107 \nexplainable ML workflow to train and test a decision tree classifier to predict endometriosis status (see 108 \nMethods) and to evaluate performance of our feature sets; 5) classifier -aligned integrated gradients (IG) 109 \nmethod to calculate gene attribution scores and 6) an LLM -powered knowledge graph (KG) approach to 110 \ninterpret biological significance of the most predictive features.  111 \n 112 \n 113 \n 114 \n 115 \nFigure 1.  End to end pipeline.  Pipeline consists of 1) identification of GEO ids for relevant 116 \nendometriosis/control datasets and LLM -powered metadata curation with human validation to ensure 117 \nconsistent and high-quality target label generation ; 2) extraction of embeddings from pre -trained frozen 118 \nRNA FMs; 3) generation of robust cohort-aware multi-fold strategy to test within vs. cross cohort learning 119 \ncapabilities of feature sets; 4) an explainable ML workflow to evaluate endometriosis classification 120 \nperformance of our feature sets; 5) classifier-aligned integrated gradients (CA-IG) method to calculate gene 121 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nattribution scores and 6) an LLM -powered knowledge graph (KG) approach to interpret biological 122 \nsignificance of the most predictive features.  123 \n 124 \nOur dataset incorporated a multi-cohort bulk RNA-seq benchmark spanning 12 independent GEO studies 125 \n(334 samples) , one of the largest collection s of independent cohorts for bulk RNA -seq endometriosis 126 \nclassification to date, and the first multi -cohort evaluation of frozen RNA foundation -model embeddings 127 \nfor this task. This approach allowed our data to be enhanced with the prior biological knowledge of the FM, 128 \nwithout the computational expense of pre -training or fine -tuning. For interpretation of the embeddings 129 \nderived from the FMs, we applied a hybrid classifier-aligned integrated gradients (CA-IG) method, which 130 \ncan produce stable, biologically plausible gene -level attributions across evaluation regimes. To identify 131 \nunexplored relationships between our most explainable genes and endometriosis in a sparsely characterized 132 \nresearch area, we applied a literature -driven knowledge graph and LLM-based inference method to 133 \nsystematically integrate and reason over published evidence, yielding a literature-grounded representation 134 \nof gene-endometriosis relationships suitable for hypothesis generation. 135 \n 136 \nThis work enabled a robust assessment of the capability of FM -based representations to generalize across 137 \nheterogenous sources, a persistent challenge situation in computational biology, with often limited 138 \navailability of large, harmonized datasets. Many studies rely on small sample sizes collected under specific 139 \nconditions, which commonly leads to individual models with good performance but poor generalizability. 140 \nThis lack of generalizability is of particular concern with understudied complex diseases like endometriosis. 141 \nBy using FM-derived embeddings in this way we aim to address this fundamental issue and understand if 142 \nrepresentations learnt from pre -trained foundation models capture biologically genuine patterns which 143 \ntranscend dataset-specific biases. 144 \n 145 \n 146 \nResults & Discussion 147 \n 148 \nNovel pipeline design enables evaluation of cross-cohort generalizability 149 \n 150 \nIn this work we developed a comprehensive end-to-end pipeline enabling inference-only use of pre-trained 151 \nRNA foundation models  for reproducible endometriosis disease prediction and gene-level interpretation 152 \nfrom bulk RNA-seq data (Fig 1). This pipeline integrates: (i ) automated download and harmonization of 153 \nGEO RNA-seq studies with LLM -assisted metadata curation and expert validation, (ii) extraction of FM 154 \nembeddings via a single forward pass through the model, (iii) cohort-aware evaluation using within-cohort 155 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nand cross-cohort strategies, (iv) consistent classical ML training and evaluation using AutoXAI4Omics 156 \ntookit, (v) explainability using SHAP for TPM models and classifier-aligned integrated gradients (CA-IG) 157 \nfor embedding models and (vi) literature -grounded knowledge-graph generation and interpretation, and 158 \npathway analysis. We assembled a multi -cohort benchmark spanning  12 independent Gene Expression 159 \nOmnibus GEO studies (334 samples ; 17,518 unified genes ) (20), by searching for human studies of 160 \nendometrial tissue related to endometriosis and restricting results to expression profiling by high -161 \nthroughput sequencing, comprising 259 endometriosis cases and 75 control s (Table 1). This multi -study 162 \nsetting captures clinically realistic heterogeneity across laboratories and protocols and enables an accurate 163 \nassessment of generalization under cohort shift.  164 \n 165 \nFrom these datasets we derived two feature representations for each sample; (i) gene-level log-transformed 166 \nTPM values, used as a baseline, and (ii) embeddings extracted from state-of-the-art ‘frozen’ RNA FMs. To 167 \nassess if FM-derived features provide an advantage over baseline features, we evaluated two classification 168 \napproaches. In approach 1, we used samples from all 12 cohorts to build the training and test sets, simulating 169 \na ‘within-cohort’ approach (see Methods). Here we expect elevated performance since we do not ask the 170 \nmodel to generalize to a completely unseen cohort, but in doing so we risk models learning patterns distinct 171 \nto each dataset . For approach 2, we designed test and train sets which consisted of distinct cohorts, 172 \nsimulating a ‘cross-cohort’ approach (see Methods). In this approach we trained the model on samples 173 \nfrom specified cohorts and tested exclusively on samples from independent cohorts which were not seen 174 \nduring training. Here we expect a drop -off in performance as the model must generalize across many 175 \ncohorts. The latter scenario reflects a common challenge in clinical and biological research, where data are 176 \ntypically distributed across many small, independent cohort. By comparing performance across these two 177 \nstrategies, we aimed to determine whether FM-derived embeddings improve robustness and generalization 178 \nin the presence of cohort heterogeneity, relative to classical gene-level features as a baseline. 179 \n 180 \n 181 \nTable 1: Summary of the cohorts used in this study. Detailing  GEO ID’s, sample counts and positive 182 \nproportions of patients or samples within each cohort with a positive diagnosis of endometriosis. 183 \n 184 \nDataset ID Dataset name (GEO name) Sample count(n)  Proportion of disease cases \n1  GSE134056              38 0.42 \n2  GSE202571              17 1 \n3  GSE205494              12 0.5 \n4  GSE135485              58 0.93 \n5  GSE212787              17 0.65 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n6  GSE135640              36 0.5 \n7  GSE226870--GSE226872   53 1 \n8  GSE157153              66 1 \n9  GSE230956              8 0.5 \n10  GSE171653              11 1 \n11  GSE190580              7 0.43 \n12  GSE243157--GSE243158   11 0 \n 185 \nWithin-cohort disease prediction  186 \n 187 \nWe first evaluated the within-cohort performance (where training and test sets include samples from the 188 \nsame cohort, though not the same patient ) for endometriosis disease classification (Table 2). To obtain a 189 \nrobust analysis of model performance variability and to allow statistical  comparison, we evaluated 190 \nperformance by performing 3-cross validation folds with 5 random-seed runs per fold per model (15 runs 191 \nin total per model), and performed hierarchical bootstrapping (10,000 iterations) with re-sampling, allowing 192 \ncalculation (for each metric) of mean, mean effect size , 95% confidence interval (CI) and two-tailed p-193 \nvalues (see Methods). 194 \n 195 \nML classification of endometriosis disease using the log transformed TPM data (baseline) resulted in a 196 \nmean weighted-F1 score of 0.86 , a class 1-F1 of 0.91 (disease), a class 0-F1 of 0.69 (no disease) and an 197 \naccuracy of 0.86 ( Table 2, Fig 2A), indicating good performance for  endometriosis classification from 198 \nmulti-cohort transcriptomic data, with more accurate prediction of disease than healthy samples (Table 2).  199 \n 200 \n 201 \nTable 2.  Within-cohort results . Predictive performance metrics for endometriosis disease prediction 202 \nusing our within-cohort strategy. Four model performance metrics are provided (weighted -F1, class0-F1, 203 \nclass1-F1 and accuracy) from classification models built using the baseline TPM features and FM -204 \nembedding features from five RNA FMs. For GeneFormer the V2 -104M-HC model was used and for 205 \nBMFM-RNA the 110m -mlm-rda-v2 model was used. For each metric, table reports the estimated mean 206 \nperformance over 10,000 bootstrap replicates and corresponding low (2.5 th) and high (97.5th) values from 207 \ncalculating the 95% confidence interval (CI) obtained via nested bootstrapping which re-sampled folds and 208 \nseed runs with replacement. All metrics are computed on the held -out test set, and ML models were built 209 \nand evaluated using AutoXAI4Omics ML pipeline (see Methods), using AdaBoost.  210 \n 211 \nModel Name Weighted F1  Class 0 – F1 Class 1 – F1 Accuracy \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n(mean [95%CI]) (mean [95%CI]) (mean [95%CI]) (mean [95%CI]) \nTPM baseline 0.86 [0.84,0.87] 0.69[0.67, 0.71] 0.91 [0.90,0.92] 0.86 [ 0.84, 0.87] \nBulkRNABert 0.89 [0.88,0.91] 0.77 [0.75, 0.79] 0.93 [0.92,0.94] 0.89 [ 0.88, 0.91] \nscFoundation 0.85 [0.81, 0.87] 0.63 [0.55,0.69] 0.91 [0.89,0.92] 0.85 [0.82, 0.88] \nscGPT 0.85 [0.81, 0.87] 0.65 [0.57,0.72] 0.91 [0.89,0.93] 0.86 [0.82, 0.89] \nGeneFormer  0.90 [0.86, 0.91] 0.76 [0.69,0.81] 0.93 [0.91, 0.94] 0.89 [0.86, 0.91] \nBMFM-RNA  0.85 [0.81,0.90] 0.66[0.55, 0.77] 0.90 [ 0.88, 0.94] 0.85 [0.70,0.81] \n 212 \n 213 \nNext, we computed the statistical significance of pairwise model performance by calculating two-tailed p-214 \nvalues from the bootstrapped samples. Among the five frozen FMs, embeddings extracted from 215 \nBulkRNABERT and Geneformer surpassed the predictive performance of the TPM baseline classifier in 216 \nall four metrics with statistical significance (P < 0.05, Supplementary Table S1A), with the weighted F1 217 \nincreasing from 0.86 for the TPM classifier to 0.89 and 0.90 for BulkRNABERT and Geneformer classifiers 218 \nrespectively (Table 2). Prediction of the healthy class (class 0) was also significantly improved using 219 \nBulkRNABERT and Geneformer embedding classifiers, increasing from a baseline-F1 of 0.69 to F1 scores 220 \nof 0.77 and 0.76 for BulkRNABERT and Geneformer classifiers respectively (Table 2, Supplementary 221 \nTable S2A, Fig 2). For the remaining three FM embedding classifiers (scFoundation, scGPT and BMFM -222 \nRNA) we saw reduced or equal model performance compared to the baseline TPM classifier for all metrics, 223 \nand for prediction of the healthy class (class 0), the scFoundation embedding classifier performed 224 \nsignificantly worse than the TPM model (P < 0.05, Supplementary Tables S2A, S3A and S4A). 225 \nSurprisingly, our top performing FM-embedding classifier was using embeddings extracted from the single-226 \ncell pretrained FM Geneformer, whereas we were using bulk RNA-seq data. This suggests that single-cell 227 \npriors transfer effectively to bulk RNA-seq.  228 \n 229 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n 230 \n 231 \nFigure 2: Within-cohort results. (a) Mean predictive performance metrics for endometriosis using our 232 \nwithin-cohort strategy for the baseline TPM features versus using frozen embeddings as from the five RNA 233 \nFMs as features. Bars show the mean across outer folds and seeds; whiskers denote 95% CIs. Metrics 234 \nreported: weighted-F1 (primary), F1 for the positive class (class 1; endometriosis), F1 for the negative class 235 \n(class 0; controls), and accuracy. (b) Effect size relative to the baseline: mean difference in score 236 \n(embedding model minus baseline) with 95% Cis. All models use identical training and evaluation within 237 \nthe series-aware within-cohort splits 238 \n 239 \nFor each FM-derived feature set classifier, mean effect-size (+ 95% CI) was also calculated which indicated 240 \nthe average difference in performance metrics from the raw baseline transcriptomic feature set model, 241 \nwhich we treated as our baseline (Fig. 2B, Supplementary Table S5). Although we observed an increase in 242 \nF1 score and positive F1 effect sizes for two FM -derived embedding classifiers, increases were moderate 243 \nin scale with very small effect sizes, and we observed a drop in model performance for the remaining three 244 \nFM embedding classifiers. Overall, within -cohort evaluation yielded high performance for both feature 245 \nrepresentation type (TPM and embeddings), with only modest and model -dependent benefits from 246 \nembeddings (Fig 2B). This pattern is consistent with a setting where cohort -specific structure is shared 247 \nbetween test and train data, and in these conditions a well-tuned classical pipeline using gene-level features 248 \nalready performs strongly, with embeddings adding limited advantage. 249 \n 250 \nFM-derived embeddings substantially improve cross-cohort disease prediction 251 \n 252 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nWe next evaluated  cross-cohort generalization, where entire cohorts are held out for testing.  For this 253 \napproach, we constructed balanced grouped cohort splits with K=3 outer folds to match the within-cohort 254 \napproach (see Methods). Individual cohorts were assigned to test folds by a constrained random search 255 \n(1,000 trials; fixed seed) that (i) enforced a minimum test size (at least 60 samples per fold) and (ii) 256 \nminimized the across-fold standard deviation of the positive fraction (disease positive samples) and of fold 257 \nsizes. We report the resultant test/train dataset assignments in Methods. Mean model performance metrics 258 \nand effect sizes were calculated as previously described (see Methods).  259 \n 260 \nWhen using our cross-cohort approach for disease prediction, the baseline classifier (log-transformed TPM) 261 \nachieved a mean weighted F1-score of 0.68, a marked decline from the F1-score of 0.86 obtained from the 262 \nwithin-cohort approach (Table 3). This loss in  performance was expected, as the trained model was less 263 \nable to generalize to samples from unseen cohorts, likely due to differences in sample collection, patient 264 \nclinical diversity, feature distribution and class-imbalances. All embedding-based classifiers also exhibited 265 \na drop in predictive performance when trained and tested using the cross-cohort vs. within-cohort approach, 266 \nhowever the degree of reduction was substantially lower in 4/5 of the embedding classifiers compared to 267 \nthe baseline classifier, with drops in mean weighted F1-score of 0.07, 0.07, 0.10 and 0.11 for GeneFormer, 268 \nscGPT, BMFM-RNA and BulkRNABert models respectively, compared to a drop of 0.18 for the baseline 269 \n(Table 3). 270 \n 271 \n 272 \nTable 3.  Predictive performance metrics from endometriosis disease prediction (cross-cohort). Four model 273 \nperformance metrics are provided (weighted-F1, class0-F1, class1-F1 and accuracy) from models using the 274 \nbaseline TPM features and FM -embedding features from 5 SOTA transcriptomic FMs. For GeneFormer 275 \nthe V2-104M-HC model was used and for BMFM-RNA the 110m-mlm-rda-v2 model was used.  For each 276 \nmetric, table reports the estimated mean performance over 10,000 bootstrap replicates and corresponding 277 \nlow (2.5th) and high (97.5th) values from calculating the 95% confidence interval (CI) obtained via nested 278 \nbootstrapping which re -sampled folds and seed runs with replacement. All metrics are computed on the 279 \nheld-out test set, and ML models were built and evaluated using AutoXAI4Omics ML pipeline (see 280 \nMethods), using AdaBoost.  281 \n 282 \nModel Name Weighted F1  \n(mean[95%CI]) \nClass 0 – F1 \n(mean [95%CI]) \nClass 1 – F1 \n(mean [95%CI]) \nAccuracy \n(mean [95%CI]) \nTPM baseline 0.68 [0.59,0.78] 0.42 [0.23, 0.53] 0.75 [0.69,0.85] 0.65 [ 0.57, 0.77] \nBulkRNABert 0.78 [0.68, 0.90] 0.53 [0.21, 0.73] 0.86 [0.82,0.91] 0.79 [ 0.71, 0.86] \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nscFoundation 0.58 [0.52, 0.63] 0.31 [0.08,0.41] 0.66 [0.55,0.79] 0.57 [0.50,0.66] \nscGPT 0.78 [0.74, 0.81] 0.50 [0.32,0.61] 0.86 [0.85,0.87] 0.79 [0.77, 0.80] \nGeneFormer  0.83 [0.77, 0.87] 0.62 [0.44,0.74] 0.89 [0.86, 0.91] 0.83 [0.78,0.87] \nBMFM-RNA 0.80 [0.71,0.86] 0.56 [0.29,0.70] 0.87 [0.83,0.91] 0.81 [0.73,0.86] \n 283 \n 284 \nIn contrast to our within -cohort approach, now when using a cross -cohort approach, all FM -embedding 285 \nclassifiers significantly outperformed the TPM -baseline model for disease classification, except for 286 \nscFoundation ( P <  0.05, Supplementary Table S1B). Geneformer and BMFM -RNA were the top two 287 \nperforming embedding models for all metrics (mean weighted F1-score, class 0 F1-score, class 1 F1-score 288 \nand mean accuracy) with mean weighted F1-scores of 0.83 and 0.80 respectively (Table 3, Fig 3, P < 0.05). 289 \nSignificant increases in performance for all embedding models apart from scFoundation compared to 290 \nbaseline were also apparent for class 1 F1 -score and accuracy metrics ( P <  0.05), while GeneFormer, 291 \nBMFM-RNA and scGPT also achieved significantly higher F1-scores for class 0 compared to baseline (P 292 \n< 0.05, Table 3, Supplementary Table S1B). For the outperforming FMs, we saw an average positive effect 293 \nsize (across all four models) compared to the baseline of 0.11 for weighted F1, class 0 F1-score and class1 294 \nF1-score and of 0.13 for accuracy (Supplementary Table S6). We saw the greatest performance 295 \nimprovement from Geneformer and BMFM -RNA embeddings, where class 0 F1 increased (compared to 296 \nTPM baseline) by an average of 0.2 and 0.14 respectively, and the weighted F1 -score increased by an 297 \naverage of 0.15 and 0.12 respectively, reflecting a substantial and significant improvement in model 298 \nperformance (Table 3, Fig 3, Supplementary Table S6). 299 \n 300 \n 301 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n 302 \n 303 \nFigure 3: Cross-cohort results (balanced grouped -cohort validation). (a) Mean performance for the 304 \n(TPM-Baseline) versus the same pipeline using frozen embeddings (BulkRNABERT, scFoundation, 305 \nscGPT, Geneformer, BMFM -MLM-RDA). Bars show the mean across outer folds and seeds; whiskers 306 \ndenote 95% CIs. Metrics: weighted -F1 (primary), F1 for the positive class (class 1), F1 for the negative 307 \nclass (class 0), and accuracy. (b) Effect size relative to the baseline: mean score difference (embedding 308 \nmodel minus baseline) with 95% CIs; the vertical dashed line at 0 marks parity (values > 0 favour 309 \nembeddings). In this setting, train and test cohorts are disjoint: each outer fold holds out a group of GEO 310 \nSeries (i.e., of cohorts) for testing, and all model selection is confined to the training side 311 \n 312 \nGene-level attribution enables interpretation of baseline TPM model 313 \n 314 \nWe investigated explainability of the within-cohort and cross-cohort baseline TPM classifiers using SHAP 315 \n(Shapley Additive Explanations) a model -agnostic attribution framework which compute s feature 316 \nimportance scores using game theory  (21). This approach allows us to identify the genes most important 317 \nfor disease prediction. Using results from all folds (n=3) and all runs (n=5) per model, for each true positive 318 \nsample (endometriosis disease positive, correctly predicted as disease) we first extracted SHAP values for 319 \nall genes and normalized across folds . To focus on disease -driving signal o ur normalization approach 320 \nincluded first setting all negative SHAP scores to zero and normalizing the remaining positive samples, 321 \nresulting in proportionate importance scores for disease-driving genes. Normalized SHAP values were then 322 \naveraged across models, runs and folds, resulting in one average SHAP value per gene per model. For the 323 \nwithin-cohort baseline TPM model (AdaBoost) across all runs and folds there were an average of  75 true 324 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \npositive samples from which the SHAP values were normalized and averaged. A total of 51 genes exhibited 325 \nnon-zero SHAP values, with the top 20 disease-associated genes shown in Table 4.  For the cross-cohort 326 \nbaseline TPM model (AdaBoost) across all runs and folds there were an average of 111 true positive 327 \nsamples from which the SHAP values were normalized and averaged with a total of 64 genes that exhibited 328 \nnon-zero SHAP values (Table 4).  329 \n 330 \n 331 \nTable 4.  Baseline TPM explainability. Top 20 positive (disease-driving) genes identified via SHAP for 332 \nTPM baseline model for the within -cohort approach and compared to the cross -cohort approach. 333 \nNormalized using a positive -only approach, then averaged across runs and folds (HGNC symbols) (Full 334 \ngene list available in Supplementary Tables S7 (within-cohort) and S8 (cross-cohort)).  335 \n 336 \nWITHIN COHORT CROSS COHORT \nGeneID score Gene name GeneID score Gene name \nEFEMP2 0.17941397 Fibulin 4-encoding gene UCN2 0.0951607 urocortin 2 \nRASSF6* 0.08950262 \nRas Association Domain \nFamily Member 6-encoding \ngene \nPCOLCE 0.09124098 \nprocollagen C-\nendopeptidase \nenhancer \nSLC39A13 0.0727763 Solute Carrier Family 39 \nMember 13 (zinc transporter) BMP1* 0.07026173 \nbone \nmorphogenetic \nprotein 1 \nTP53I13 0.064122 Tumor Protein P53 Inducible \nProtein 13 LLGL1* 0.05474552 \n LLGL scribble \ncell polarity \ncomplex \ncomponent 1 \nBMP1* 0.05455228 Bone Morphogenetic Protein 1-\nencoding gene MAL 0.05190146 \nmal, T cell \ndifferentiation \nprotein (MAL \nblood group) \nLCN2 0.04278986 Lipocalin 2 -encoding gene CYP8B1 0.05036643 \ncytochrome \nP450 family 8 \nsubfamily B \nmember 1 \nSTEAP4 0.03793841 STEAP4 Metalloreductase RASSF6* 0.04911975 \nRas association \ndomain family \nmember 6 \nVPS16 0.03228817 Vacuolar Protein Sorting \nProtein 16-encoding gene TSPAN4 0.04771203 tetraspanin 4 \nCDC37 0.0305804 Cell Division Cycle 37, HSP90 \nCochaperone AGRN 0.03964136 agrin \nSAMD11* 0.03008721 Cell division cycle 37 homolog SAMD11* 0.03699761 \nCell division \ncycle 37 \nhomolog \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nMXRA8* 0.02810625 Matrix remodeling-associated \nprotein 8 ARFGAP1 0.03572216 \nARF GTPase \nactivating \nprotein 1 \nOAF 0.02796728 Out at first homolog IDO1 0.03315175 \nindoleamine \n2,3-\ndioxygenase 1 \nPDK4 0.02608833 pyruvate dehydrogenase kinase \n4 \nCCDC141 0.03049321 \ncoiled-coil \ndomain \ncontaining 141 \nSURF4 0.02398875 surfeit 4 CD276 0.02285808 CD276 \nmolecule \nHPS1 0.02067466 HPS1 biogenesis of lysosomal \norganelles complex 3 subunit 1 MLXIPL 0.02278676 \nMLX \ninteracting \nprotein like \nTMEM273 0.01807645 transmembrane protein 273 MXRA8* 0.02076158 \nmatrix \nremodeling \nassociated 8 \nTAX1BP3 0.01622053 Tax1 binding protein 3 ITGAM 0.01675609 integrin subunit \nalpha M \nLLGL1* 0.01517964 LLGL scribble cell polarity \ncomplex component 1 \nCHST14 0.01642441 \ncarbohydrate \nsulfotransferase \n14 \nADAMTS2 0.013455 ADAM metallopeptidase with \nthrombospondin type 1 motif 2 \nTRAF7 0.01502836 \nTNF receptor \nassociated \nfactor 7 \nLYZ 0.01336511 lysozyme TMEM129 0.01500533 \ntransmembrane \nprotein 129, E3 \nubiquitin ligase \n 337 \n 338 \nWe observed little overlap between the TPM-baseline explainability results for the within-cohort vs cross-339 \ncohort approaches, with only 5/20 of the top 20 explainable genes in the within-cohort approach also within 340 \nthe top 20 for the cross-cohort approach (Table 4). The top four most important genes for the within-cohort 341 \nand cross-cohort approaches were completely distinct, being derived from different biological pathways 342 \n(WikiPathways), indicating that different biological signals were being captured by the within and cross-343 \ncohort approaches (Supplementary Figure 1, Table 4).  This variability in the predictive features between 344 \nmodels is concerning for researchers using interpretation of ML models to guide biomarker or target 345 \ndiscovery for disease and could partially explain our drop in performance for cross-cohort analysis if the 346 \npredictive features from the within-cohort performance do not generalize in the cross-cohort context i.e., 347 \ncohort specific predictors are being learned in the within-cohort context.  348 \n 349 \nThe five overlapping genes between the within and cross-cohort approaches were RASSF6, BMP1, MXRA8, 350 \nSAMD11 and LLGL1. RASSF6 encodes a Ras-effector protein which binds activated Ras and is involved in 351 \ninduction of apoptosis . It’s a member of a tumor repressor network which is highly conserved , and it’s 352 \ninactivation or suppression has been linked to several diseases  including cancer  (22). BMP1 (bone 353 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nmorphogenic protein 1) is a metalloprotease which regulates formation of extracellular matrix (ECM) by 354 \ncleaving/processing of various precursor proteins. Mutations in this gene have been associated with 355 \nskeletal/connective tissue disorders (23). The role of MXRA8 (Matrix Remodeling Associated 8- encoding 356 \ngene) is less established, but it is thought to be involved in glial blood -brain barrier functioning and has 357 \nbeen associated with glioma progression and a suggested therapeutic target for glioma  (24). MXRA8 is 358 \nknown to play roles in inflammation and cell signaling that could be potentially relevant to endometriosis. 359 \nLong non-coding RNA transcribed from the SAMD11 region has been functionally linked to endometrial 360 \nbiology (decidualization) (25). LLGL1 is a well -known tumor suppressor gene involved in maintaining 361 \nepithelial cell polarity and tissue architecture in various cancers, including endometrial cancer (26).    362 \n 363 \nClassifier-aligned integrated gradients (CA -IG) allows FM embedding  gene–level 364 \ninterpretation 365 \n 366 \nNext, we sought to understand if and how model explainability differed when using FM embeddings for 367 \ndisease classification. Because the downstream classifier is not jointly fine-tuned with the encoder, standard 368 \ngradient attribution cannot be directly applied. We therefore developed classifier -aligned integrated 369 \ngradients (CA -IG), an attribution approach that assigns importance to FM gene tokens by computing 370 \nIntegrated Gradients (IG) through the frozen encoder along a classifier-relevant direction in the embedding 371 \nspace. This direction is defined locally for each sample using SHAP scores from the downstream classifier 372 \nwith respect to the FM embedding, ensuring that gene -level attributions reflect the features used by the 373 \nclassifier. Importantly, SHAP is used solely to define this classified -aligned weighting of embedding 374 \ndimensions, while IG performs the attribution itself. This formulation avoids backpropagation through the 375 \nnon-differentiable classifier, substantially reducing computational cost while preserving biological 376 \ninterpretability by providing gene-level explanations (see Methods).  377 \n 378 \n 379 \nTable 5. FM approach. Top 20 explainable genes identified via CA-IG for Geneformer embedding model 380 \nfor the within-cohort and cross-cohort approach. SHAP scores normalized using a positive-only approach, 381 \nthen averaged across runs and folds (HGNC symbols). Asterisks (*) represent  genes found  by both 382 \napproaches. (Full gene list available in Supplementary Tables S9 (within-cohort) and S10 (cross-cohort)).  383 \n 384 \nWITHIN COHORT CROSS-COHORT \nGeneID CA-IG score Gene name GeneID CA-IG score Gene name \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nDDIT3* 0.01254973 DNA Damage Inducible \nTranscript 3 DDIT3* 0.00966168 DNA Damage \nInducible Transcript 3 \nLRRC3C* 0.00604412 leucine rich repeat \ncontaining 3C \nLRRC3C* 0.0058623 leucine rich repeat \ncontaining 3C \nTBC1D3F* 0.00598382 TBC1 domain family \nmember 3F \nTBC1D3F* 0.00584974 TBC1 domain family \nmember 3F \nOR1J2* 0.00528166 olfactory receptor family 1 \nsubfamily J member 2 \nOR1J2* 0.00370744 \nolfactory receptor \nfamily 1 subfamily J \nmember 2 \nFRG2* 0.00337272 FSHD region gene 2 FRG2* 0.00330829 FSHD region gene 2 \nTBC1D3E* 0.00240315 TBC1 domain family \nmember 3E \nMEPCE* 0.00101443 methylphosphate \ncapping enzyme \nPPIAL4F* 0.0015164 peptidylprolyl isomerase A \nlike 4F \nTBC1D3E* 0.00099277 TBC1 domain family \nmember 3E \nMEPCE* 0.00144943 methylphosphate capping \nenzyme \nTBC1D3K* 0.0009618 TBC1 domain family \nmember 3K \nOR51B5 0.00127887 olfactory receptor family \n51 subfamily B member 5 \nTBC1D3* 0.00094298 TBC1 domain family \nmember 3 \nTBC1D3K* 0.00124146 TBC1 domain family \nmember 3K \nRSC1A1* 0.00092813 regulator of solute \ncarriers 1 \nTBC1D3C* 0.00122356 TBC1 domain family \nmember 3C \nPPIAL4F* 0.00089612 peptidylprolyl \nisomerase A like 4F \nRABGAP1* 0.00120914 RAB GTPase activating \nprotein 1 \nRABGAP1* 0.00089356 RAB GTPase \nactivating protein 1 \nRSC1A1* 0.00119487 regulator of solute carriers \n1 \nDERPC 0.00088404 \nDERPC proline and \nglycine rich nuclear \nprotein \nTBC1D3* 0.00118897 TBC1 domain family \nmember 3 \nSPDYE14* 0.00075791 \nspeedy/RINGO cell \ncycle regulator family \nmember E14 \nTBC1D3H* 0.00116252 TBC1 domain family \nmember 3H \nTBC1D3C* 0.00074984 TBC1 domain family \nmember 3C \nSPDYE10* 0.00109991 \nspeedy/RINGO cell cycle \nregulator family member \nE10 \nPCDHGC5* 0.00073012 protocadherin gamma \nsubfamily C, 5 \nSPDYE15 0.00100317 \nspeedy/RINGO cell cycle \nregulator family member \nE15 \nTBC1D3H* 0.0007123 TBC1 domain family \nmember 3H \nPCDHGC5* 0.00097082 protocadherin gamma \nsubfamily C, 5 \nPIGY 0.00069504 \nphosphatidylinositol \nglycan anchor \nbiosynthesis class Y \nPARP2* 0.00096208 poly(ADP-ribose) \npolymerase 2 \nPARP2* 0.0006938 poly(ADP-ribose) \npolymerase 2 \nSPDYE14* 0.00093309 \n \nspeedy/RINGO cell cycle \nregulator family member \nE14 \nSPDYE10* 0.00068962 \n \nspeedy/RINGO cell \ncycle regulator family \nmember E10 \n 385 \n 386 \n 387 \nUsing our CA -IG approach to extract explainability for the FM -embedding classifiers, we observed a 388 \nstriking conservation of disease -associated genes between the within -cohort and cross -cohort settings 389 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n(Table 5).  Of the top 20 genes identified in the within -cohort analysis, 18 overlapped within the cross -390 \ncohort top 20, and conversely 18/20 of the top cross-cohort genes were also present in the within-cohort list 391 \n(Table 5). Notably, the top five genes were identical and ranked in the same order for both strategies, 392 \nindicating highly consistent relative importance across evaluation regimes (Fig 4). This contrasts sharply 393 \nwith the TPM-based models, where predictive features varied substantially between within and cross cohort 394 \nanalyses. The stability of the embedding -derived attributions parallels the reduced performance drop 395 \nobserved for the FM-based classifiers under cohort shift, suggesting that the pretrained representations may 396 \nprioritise transferable disease-related signal over cohort-specific effects.  397 \n 398 \nThe five most predictive genes across both settings were DDIT3 (DNA Damage-inducible transcript 3 -399 \nencoding gene), LRRC3C (Leucine Rich Repeat Containing 3C-encoding gene), TBC1D3F (TBC1 Domain 400 \nFamily Member 3F), OR1J2 (Olfactory Receptor Family 1 Subfamily J Member 2 -encoding protein) and 401 \nFRG2 (FSHD Region Gene 2-encoding protein) (Fig 4). 402 \n 403 \n 404 \n 405 \n 406 \nFigure 4: Conserved gene-level attributions from Geneformer embeddings. Ten highest–scoring genes 407 \nfor the frozen Geneformer  encoder, computed with Integrated Gradients (PAD baseline; Methods) and 408 \naggregated over true-positive test samples, runs/seeds, and outer folds. Panels show the within-cohort (left) 409 \nand cross-cohort (right) strategies. Values are per -sample positive-only, L1-normalised attribution scores 410 \n(positive SHAP values kept only), with larger indicates more disease-driving evidence. 411 \nWITHIN-COHORT CROSS-COHORT     \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n 412 \nThe most important gene for endometriosis disease prediction using both the within-cohort and cross-cohort 413 \napproaches was DDIT3, which encodes a transcription factor which  mediates response to cellular stress , 414 \nwith a main role in promoting apoptosis during conditions of severe stress. Endoplasmic reticulum (ER) 415 \nstress is being increasingly recognized as playing a major role in the pathogenesis of endometriosis (27,28), 416 \nand DDIT3/CHOP is a crucial transcription factor activated by ER stress which has also been shown to 417 \nfurther upregulate ER stress and apoptosis in endometriosis disease (29). DDIT3 promotes apoptosis (30), 418 \nwhich may drive the degeneration of ectopic tissue.  DDIT3 has been strongly linked to a variety of other 419 \ndiseases including inflammatory diseases, cancer, neurodegenerative diseases and metabolic diseases  420 \n(31,32).  Its’ overexpression has been correlated with poor breast cancer prognosis and implicated in disease 421 \nprogression via its role in cellular stress response pathways (33). DDIT3 plays key roles in inflammatory 422 \ndiseases by activation of several inflammatory pathways including NF-kb pathway (34) and is upregulated 423 \nby cytokines TNF -a and IL-1B. DDIT3 overexpression has also been linked to reduced antiviral innate 424 \nimmune activity, encouraging viral replication and disease (35). This gene is a strong candidate for further 425 \nexploration and validation.  426 \n 427 \nThe second ranked gene, LRRC3C, encodes a (predicted) membrane -located protein containing leucine -428 \nrich repeats motifs (LRR), important for protein-protein interactions. Although mechanistic studies linking 429 \nLRRC3C with disease remain limited, LRRC3C has been implicated as important in inflammatory disease, 430 \nparticularly in inflammatory bowel disease susceptibility via its impact on cell proliferation and apoptosis 431 \n(36).  432 \n 433 \nThe protein encoded by the third most important gene TBC1D3F is one of eight paralogs (A-H) belonging 434 \nthe TBC1D3 family. Members of this family of molecules have been found to be over expressed in several 435 \ncancer types (prostate, breast, pancreatic) (37) and have been associated wi th poor prognosis of kidney 436 \nrenal cell carcinoma (38). TBC1D3F was found to be overexpressed in the epithelial ovarian cancer cell 437 \nlines and primary tumors compared to normal tissues (39). Alongside TBC1D3F, several members of the 438 \nTBC1D3 family are present in the top 20 list of most important genes for endometriosis disease prediction, 439 \nincluding TBC1D3E, TBC1D3K and TBC1D3RC. The TBC1D3 locus is frequently associated with cancer, 440 \namplified in 15% of prostate cancers , and the TBC1D3 cluster is extremely structurally variable (40,41). 441 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nThe amplified presence of genes from the TBC1D3 family ( 6/20 of the topmost important genes for both 442 \nwithin and cross-cohort models) makes this family a strong candidate for experimental validation.   443 \n 444 \nOlfactory receptors (ORs) are expressed in several non -olfactory tissues including the endometrium and 445 \nare linked to inflammation through specific signaling pathways. Although information on OR1J2 is 446 \nrelatively sparse, other ORs, for example, OR5H2, have been identified as a potential target in endometrial 447 \ncancer biology related to the IGF1 signaling pathway, and may be involved in cancer protection in some 448 \nindividuals (42).  449 \n 450 \nLLM-driven knowledge graph inference reveals enriched biological processes linked to 451 \nendometriosis 452 \n 453 \nHaving first investigated the individual functional roles of the top genes and their potential relevance 454 \nendometriosis disease pathology, we next sought to characterise  the shared molecular landscape of these 455 \ngenes to determine whether they shared any links or associations with each other and/or with endometriosis. 456 \nBecause the molecular mechanisms of endometriosis remain poorly characterised, we employed a 457 \nliterature-derived knowledge-graph and LLM-inference approach to do this, which may reveal relationships 458 \nwhich traditional bioinformatic gene-centric analysis may miss. For each of the five most important genes 459 \nfor endometriosis prediction (DDIT3, LRRC3C, TBC1D3F, OR1J2 and FRG2), we systematically retrieved 460 \nall PubMed-indexed articles in which it was mentioned and identified semantic relationships between that 461 \ngene and identified biomedical entities such as disease, gene, pathway, proteins etc. All entities and 462 \nrelationships were represented as entity -entity-relation triples and integrated into a unified bipartite 463 \nknowledge graph.    464 \n 465 \nThe resulting knowledge graph contained 312,598 triples linking 81,219 distinct entities of which 9,835 466 \nwere proteins, 8,669 were diseases, 5,800 were genes, 4,698 were biological process and 1,367 were drugs 467 \n(Full list entity types in Supplementary Table S11). The graph contains 33,092 distinct relation types with 468 \nthe most common being co-expression (4,910 triples) and regulation (4,095 triples) (Full list relation types 469 \nin Supplementary Table S12). First, we investigated the similarity of connected entities and relations 470 \nbetween the five important genes by calculating i) normalized shared -entity overlap, and ii) relation-style 471 \nsimilarity for each gene pair (see Methods). As expected, most gene pairs exhibited relatively low 472 \nnormalized shared-entity overlap, suggesting genes tended to engage with largely non -overlapping entity 473 \nsets (Supplementary Table S13). Due to the vastly larger number of triples extracted for gene DDIT3 than 474 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nthe other important genes, 27%-43% of unique entities from genes LRRC3C, TBC1D3F, OR1J2 and FRG2 475 \nwere shared with DDIT3 (Supplementary Table S13).  476 \n 477 \nTo investigate how our five important genes may converge on common biological mechanisms linked to 478 \nendometriosis disease pathology, we reduced the multi-gene knowledge graph to a compact JSON structure 479 \ncontaining shared entities, shared relations and inter -gene similarity score, and provided as structured 480 \ncontext to an LLM to suggest biologically plausible mechanisms supported by the evidence. Although the 481 \nfive genes have distinct primary functions, the LLM results highlighted their roles in several convergent 482 \nbiological themes linked to endometriosis. ER stress is being increasingly recognized as playing a major 483 \nrole in the pathogenesis of endometriosis (27,28), and is enhanced by DDIT3/CHOP in endometriosis (29). 484 \nThe LLM results correctly characterize DDIT3s potential role in endometriosis disease and links it to the 485 \nrole of TBC1D3F in vesicle-mediated cytokine release (43), suggesting a shared influence on inflammatory 486 \nand stress -response pathways in endometriosis. LLM results also suggest a role of TBC1D3F in 487 \nendometriosis pathogenesis via its role in dysregulated vesicle trafficking, which can alter secretion of 488 \nproinflammatory cytokines and growth factors that promote ectopic endometrial implantation and 489 \ncontribute towards the pro -inflammatory milieu. This is supported by the literature which indicate that 490 \nextracellular vesicles play a role in endometriosis disease pathogenesis (44). LLM results suggest a joint 491 \ncontribution of genes LLRC3C and FRG2 to endometriosis pathogenesis through processes related to 492 \nenhanced adhesion and fibrosis, mediated by extracellular -matrix remodeling ( LLRC3C) and fibroblast 493 \nproliferation/fibrosis (FRG2).  In the literature LRRC3C is only partially characterized with predicted roles, 494 \nhowever LRR (leucine-rich repeat) proteins have been implicated in cell-adhesion and extra-cellular-matrix 495 \n(ECM) remodeling (45,46). Enhanced adhesion of endometrial stomal cells to peritoneal surfaces and 496 \naltered ECM are key steps in endometrial lesion establishment (47,48).   497 \n 498 \nBMFM attribution rankings  for GSEA  reveals enriched gene pathways linked to 499 \nendometriosis pathophysiology 500 \n 501 \nTo further investigate associations between the model -derived attribution scores and potential 502 \nendometriosis disease mechanisms, we performed gene set enrichment analysis (GSEA) from the cross -503 \ncohort embedding descending gene list ranked by CA -IG score (see Methods). This analysis identified 504 \nseveral positively enriched pathways (FDR < 0.05), indicating a significant association with disease 505 \nclassification. Enriched disease pathways were predominantly related to inflammatory cytokine signaling, 506 \ncell stress and cell survival, mechanisms consistent with the current understanding of endometriosis disease 507 \npathophysiology (Supplementary Figure 1, Supplementary Tables S14 -S17). Within cytokine signaling, 508 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \npathways responsible for positive regulation of type II interferon (IFN-gamma), interleukin-17 (IL-17) and 509 \ninterleukin-8 (IL-8) were all significantly enriched amongst top ranking genes (FDR < 0.05, Supplementary 510 \nFigure 1). Recent evidence has suggested IL -17 may be an important mediator of endometriosis disease, 511 \nwith elevated levels observed across several endometriosis studies (49,50) and sample types including 512 \nendometriotic lesions (51) and peritoneal fluid (52) from women with endometriosis.  Furthermore, higher 513 \nlevels of IL -17-producing cells have been found in the peripheral blood of women with endometriosis, 514 \naccompanied by increased IL -17 protein expression in their endometriotic tissue (53), indicating both 515 \nsystemic and local involvement of IL -17 in endometriosis disease. IL -8 is also frequently elevated in the 516 \nperitoneal fluid of women with endometriosis, and in vitro studies have demonstrated IL-17 stimulation of 517 \nIL-8 secretion from endometriotic stomal cells, linking both pro -inflammatory cytokines to sustenance of 518 \nthe pro-inflammatory milieu of endometriosis (54). Elevated levels of type II interferon (IFN-gamma) have 519 \nalso been detected in endometriotic tissue and peritoneal fluid of women with endometriosis, supporting a 520 \nrole for IFN -gamma mediated inflammatory and immune dysregulation in endometriosis disease 521 \npathophysiology (55,56), however some studies have reported conflicting results (57). Alongside cytokine 522 \nsignaling, GSEA indicated several enriched pathways associated with disease were linked to cell stress and 523 \ncell survival mechanisms, including endoplasmic reticulum (ER) stress, unfolded protein response (UPR) 524 \nand regulation of apoptosis signaling pathways, highlighting altered cell stress and survival programs in 525 \nendometriosis (Supplementary Figure 1). ER stress and UPR activation have been shown in endometriotic 526 \ntissue, and is thought to reflect a response of ectopic cells to the pro -inflammatory milieu (28,58,59). 527 \nEnrichment of pathways regulating ER stress-induced apoptosis is consistent with a dysregulated balance 528 \nbetween stress responses and programmed cell death, a process which has been widely implicated in 529 \nendometriotic lesion persistence via survival of endometrial cells via evasion of apoptosis (60,61). Together 530 \nwith existing literature, these findings support the interpretation that stress-altered regulation of cell survival 531 \npathways represent a prominent feature of endometriosis pathophysiology.  532 \n 533 \nConclusions  534 \n 535 \nIn this study, we leveraged a multi-cohort bulk RNA -seq dataset to systematically assess the benefit of 536 \nRNA foundation model (FM) embeddings as feature representations for endometriosis disease prediction 537 \nusing a downstream classifier.  We evaluated the benefits of these FM-derived features using a cohort aware 538 \nstrategy and found that FM-derived representations consistently and substantially outperformed gene-level 539 \nTPM features under a cross -cohort evaluation  approach, where models were required to generalize to 540 \nentirely unseen datasets. This setting reflects a common challenge in clinical transcriptomics and biomarker 541 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \ndiscovery, and our results demonstrate that pre-trained RNA representations improve robustness to cohort 542 \nheterogeneity.  543 \n 544 \nA central contribution of this work is the introduction of classifier-aligned integrated gradients (CA-IG), a 545 \ncomputationally efficient approach for interpreting gene-derived FM embeddings. While IG benefits from 546 \na GPU for speed, it substantially lowers the computational burden for interpretability of FMs and makes 547 \nanalyses feasible with modest resources. Applying CA -IG to Geneformer embeddings revealed highly 548 \nstable gene-level attributions across both cohort-aware strategies, in contrast to TPM-based models where 549 \npredictive genes varied substantially between cross and within cohort approaches. Embedding models 550 \nconsistently prioritised a conserved set of genes, supporting the notion that pretrained encoders can capture 551 \ngeneralisable biological knowledge rather than cohort -specific noise. Pathway -level analyses further 552 \nrevealed convergence on inflammatory cytokine signaling, cell stress and cell survival mechanisms, 553 \nprocesses strongly implicated in endometriosis disease pathophysiology.  554 \n 555 \nSome limitations remain within this work. First is the small size of several of the individual cohorts, which 556 \nincreased variance in train/test splits. We addressed this using dataset-aware splits, a grouped-cohort outer 557 \nCV, a nested inner CV for model selection, and multiple seeds with bootstrap confidence intervals.  Second 558 \nwas the imbalanced nature of the full cohort, with comparatively few controls. This imbalance likely 559 \ncontributed to the reduced predictive performance of class 0 compared to class 1 and is a common issue in 560 \ndisease/control classification . Our third limitation was associated with inconsistent reporting of key 561 \ncovariates due to small and varied sample sizes (e.g., menstrual -cycle phase, lesion subtype), which 562 \nprevented adjustment for potential confounders.  563 \n 564 \nRNA FMs offer an accessible and powerful solution for improving the robustness of disease classification 565 \nin transcriptomics. Here we demonstrated this for endometriosis, a condition hampered by diagnostic delays 566 \nand biological heterogeneity, where this CPU-friendly approach provided more generalizable predictions 567 \nand stable, interpretable biomarkers, laying a practical foundation for future application to other diseases 568 \nand use-cases, leveraging large-scale biological priors in clinically focused research.  569 \n 570 \n 571 \nMethods 572 \n 573 \nDatasets and pre-processing 574 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n 575 \nWe assembled a multi -cohort bulk RNA-seq corpus from the Gene Expression Omnibus (GEO) (20) by 576 \nsearching for human studies of endometrial tissue related to endometriosis and restricting results to 577 \nexpression profiling by high-throughput sequencing. Single-cell, spatial transcriptomic studies, and studies 578 \nwith fewer than 10 samples (pre -filtering) were excluded. We used the GEO identifiers for our selected 579 \nstudies to download the associated transcriptomic data for these studies and firstly, to merge and normalise 580 \nthe gene expression count data (TPM) into a unified data matrix, and then to extract harmonised metadata 581 \nfor the same samples i.e., in this case we used an LLM to yield  a label for each sample of endometriosis 582 \ndisease yes (class 1) or control/no (class 0). Labels were derived from the harmonised metadata and then 583 \nfurther checked and curated by an expert molecular biologist to remove ambiguous cases (e.g., endometrial 584 \ncancer, unclear disease status). After filtering, the final corpus comprised 12 studies and 334 samples (259 585 \nendometriosis, 75 healthy controls). GEO Series IDs were retained as a grouping variable for subsequent 586 \ncohort-aware cross-validation. Raw counts were converted to transcripts per million (TPM) using GRCh38 587 \n(p13) gene lengths and log-transformed (log(1+TPM)).  Gene identifiers were harmonised across cohorts, 588 \nand all matrices were aligned to a common 17,518-gene universe. 589 \n 590 \nFold construction 591 \n 592 \nWithin-cohort  593 \n 594 \nFor the within -cohort analysis, we used a dataset -aware, stratified 3 -fold split approach. To ensure each 595 \ndataset was represented in the training set for each fold, we first selected one sample per GEO series or 596 \ndataset as the ‘anchor’ sample and pinned all anchors to the training set in every fold. We then applied 597 \nstratified-K-fold (K=3, shuffled, fixed seed) to the non -anchor samples only, preserving the global class 598 \nratio in each test fold. In fold k the training set was the union of all anchors and the non-anchor samples not 599 \nassigned to the k-th test split; the test set was the corresponding stratified k-th split. This procedure yields 600 \nzero sample overlap between train and test, keeps test-set class balance approximately intact, and enforces 601 \nan in-domain condition in which any series or dataset appearing in test also has at least one sample (its 602 \nanchor) in train. Fold assignments and per-fold counts are reported in Table 6. 603 \n 604 \nCross-cohort  605 \n 606 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nFor the cross-cohort approach, we constructed balanced grouped-cohort splits with K=3 outer folds to match 607 \nthe within-cohort approach. Entire datasets (GEO Series IDs) were assigned to test folds by a constrained 608 \nrandom search (1,000 trials; fixed seed) that (i) enforced a minimum test size (at least 60 samples per fold) 609 \nand (ii) minimised the across-fold standard deviation of the positive fraction and of fold sizes. This resulted 610 \nin test/train dataset assignments as reported in Table 6. 611 \n 612 \n 613 \nTable 6. Fold construction for the within-cohort and cross-cohort folds.  614 \n 615 \nWithin-\ncohort \nFolds \n Train (n)    Test (n)   \nPositive \nproportion \n(train)   \nPositive proportion \n(test)   Dataset IDs  \n                  \n1  226  108  0.78  0.77  1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12  \n2  227  107  0.78  0.77  1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12  \n3  227  107  0.78  0.78  1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12  \nCross-\ncohort \nFolds \n Train (n)    Test (n)   \nPositive \nproportion \n(train)   \nPositive proportion \n(test)   Dataset IDs  \n      \n1  237  97  0.78  0.77  6,7,9  \n2  228  106  0.77  0.78  5,8,12  \n3  203  131  0.78  0.77  2,4,10,11  \n 616 \n 617 \nExtraction of Foundation Model Embeddings 618 \n 619 \nWe evaluated five pretrained encoders spanning bulk and single -cell transcriptomic pretraining: 620 \nBulkRNABERT (18), scFoundation (17),  scGPT (16),  Geneformer (15)  and BMFM-MLM-RDA (19). 621 \nThese models learn transferable structure from large expression corpora and produce compact, context -622 \naware representations (embeddings) for downstream prediction.  623 \n 624 \nFor each sample we performed feature harmonisation as follows: (i) mapping the TPM vector to the 625 \nencoder’s gene vocabulary (e.g. HGNC to Ensembl); (ii) dropping genes absent from that vocabulary; (iii) 626 \nconsolidating duplicate symbols by averaging their TPM values; and (iv) constructing model -specific 627 \nfeature sets (tokens). To extract embeddings per sample, we then executed a single forward pass with the 628 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nencoder in evaluation mode (frozen weights) and obtained a fixed -dimensional embedding per sample. 629 \nEmbeddings were computed once per sample and cached before any model selection; identical splits and 630 \nseeds were reused across models for paired comparisons. We intentionally avoided additional downstream 631 \nscaling or feature selection on embeddings to preserve the pretrained geometry. When multiple HGNC 632 \nsymbols mapped to the same Ensembl identifier, we averaged their TPM values by grouping them by their 633 \nshared Ensembl identifier prior to input construction. We chose this approach as averaging preserves scale 634 \nand avoids artificial inflation for genes with many aliases, which could otherwise bias downstream models. 635 \nModel-specific approaches are detailed below: 636 \n 637 \nGeneformer  638 \nGeneformer is a transformer pretrained on 104M single-cell transcriptomes with rank–value encodings and 639 \nmasked-token objectives  (15). We mapped TPM(log+1) matrices to Ensembl IDs, tokenised  with the 640 \nofficial toolkit, and extracted the penultimate layer. Mean pooling excluded special tokens (CLS/EOS, 641 \nPAD) producing a 768-dimension embedding per sample.  642 \n 643 \nBulkRNABERT  644 \nThis encoder-only model was pretrained on bulk RNA -seq (GTEx/ENCODE/TCGA) with a masked-LM 645 \nobjective over discrete expression tokens and permutation-invariant gene order (18). TPM(log+1) matrices 646 \nwere harmonised to the accepted Ensembl set and ordered to the model list. We extracted embeddings from 647 \nthe final layer mean-pooled across tokens to yield a 256-dimensional embedding per sample. 648 \n 649 \nscFoundation 650 \nScFoundation is an asymmetric encoder–decoder pretrained on >50 million human single-cell profiles (62). 651 \nTPM(log+1) vectors were aligned to the 19,264 -gene vocabulary (missing genes padded with zeros). The 652 \nofficial routine returns a 3,072 -dimensional vector comprising four 768 -dimensional pools (1: last token, 653 \n2: penultimate token, 3: max pooled, 4: mean pooled). We split this file and used the last 768 dimensions 654 \nfor all further analysis (mean-pooled 768-dim block). 655 \n 656 \nscGPT 657 \nScGPT is a generative transformer pretrained on 33 million single -cell profiles with masked -attention 658 \nobjectives, jointly learning gene and cell embeddings (16). TPM(log+1)  matrices were aligned to the model 659 \nvocabulary and embedded, yielding a 512-dim vector per sample. 660 \n  661 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nBMFM-MLM-RDA 662 \nWe used the MLM -RDA variant of BMFM: a BERT -base encoder, pretrained with read -depth–aware 663 \n(RDA) down sampling pretrained via expression masking and RDA downsampling 664 \n(https://huggingface.co/ibm-research/biomed.rna.bert.110m.mlm.rda.v1, (19)). We used TPM (log+1) data 665 \nfor a single forward pass through the frozen embedder; the final -layer [CLS] vector (768 -dimensional) 666 \nserved as the sample embedding. 667 \n 668 \nExplainable Machine Learning workflow 669 \n 670 \nTrain and test protocol  671 \nWe used the open-source explainable-AI tool AutoXAI4Omics (https://github.com/IBM/AutoXAI4Omics) 672 \nto run a robust classical explainable ML analysis on both TPM and embedding inputs for both the within -673 \ncohort and cross-cohort strategies, ensuring consistency, transparency, and reproducibility across model 674 \ncomparisons (63). We chose the AdaBoostClassifier model as our downstream light-weight classifier (64).   675 \n 676 \nWe evaluated two input data types within the same AutoXAI4Omics classical -ML pipeline for 677 \nendometriosis classification: (i) log TPM values and (ii) foundation-model (FM) embeddings derived from 678 \nlog TPM values. We adopted two strategies; i) within-cohort, where train and test contain samples from the 679 \nsame GEO dataset ID, and ii) cross -cohort, where the held -out test fold comprised unseen datasets not 680 \nrepresented within the training set. The outer loop used K=3 folds; in each outer fold, models were trained 681 \non K-1 folds and evaluated once on the held-out fold. Within the fold's training set, hyperparameters were 682 \ntuned by stratified 5-fold cross-validation, selecting the configuration that maximised weighted-F1 (primary 683 \nmetric), then refitting that configuration on the full training data before testing. To assess robustness, each 684 \nouter split was repeated with five different random seeds; we report means across folds and seeds. All 685 \npreprocessing (scaling, feature transforms, encoders and feature selection) was fitted strictly on the training 686 \nside and applied unchanged to validation/test to avoid leakage.  687 \n 688 \nTPM data preparation for ML 689 \nFor all AutoXAI4Omics runs with TPM data we performed feature selection using the following two inbuilt 690 \nmethods; (i) first a variance filter removed near -constant genes (threshold = 1); and (ii) secondly a 691 \nunivariate ANOVA test (F-statistics) was used to select the top k genes, where k was chosen automatically 692 \nfrom a base-10 grid with a minimum of 25 features (no explicit maximum). For each candidate k, a small 693 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nrandom forest was trained on the selected genes and evaluated by F1 score. From the set of values yielding 694 \nhigh performance, we selected the smallest k, prioritising a smaller feature set while retaining performance.  695 \n 696 \nFM-embedding data preparation. 697 \nFor all FM vectors we disabled standardisation and feature selection, treating embeddings as fixed 698 \nrepresentations to avoid distorting pretrained geometry. 699 \n 700 \nPost-AutoXAI4Omics Bootstrapping and significance calculation  701 \n 702 \nTo obtain a robust and sound estimate of all models performance variability and allow statistically valid 703 \ncomparison of models, we calculated confidence intervals for fold -averaged metrics via hierarchical 704 \nbootstrapping, with 10,000 iterations. For each model we had performance metrics from 3 folds and 5 seed 705 \nruns per fold, totaling 15 runs per model. Each bootstrap iteration consisted of drawing 3 folds with 706 \nreplacement, and per selected folds drawing 5 seed runs with replacement. For each metric ‘fold means’ 707 \nwere calculated which were then used to calculate ‘overall means’ across all folds. The same approach was 708 \nused to calculate the effect size (while preserving paired seeds across models), defined as the difference in 709 \nmodel metrics to the baseline model (TPM input, not embeddings). Bootstrap estimates of the mean and 710 \neffect size over 10,000 iterations were used to construct 95% confidence intervals. To test whether 711 \nperformance of paired models differed significantly we leveraged the same nested -bootstrap samples that 712 \nwe used for confidence intervals. In each of the 10,000 bootstrap iterations we paired the corresponding 713 \nruns of two models (same index), computed the difference in performance metric for each fold, and 714 \naveraged these differences over the 3 folds. This yielded one bootstrap difference value per metric per 715 \niteration. Under the null hypothesis that the models have equal performance, the distribution of differences 716 \nwould be centred on zero (no difference). We therefore we calculate two -tailed p-values as follows:  i) 717 \ncompute the proportion of bootstrap differences which lie above zero, and the proportion that lie below 718 \nzero, ii) take the smaller proportion (representing the least likely direction of difference) and iii) we double 719 \nthis to get a 2 -tailed p-value. This approach was applied to all pairwise comparisons (baseline/vanilla vs. 720 \nall BMFMs) for the metrics F1-score, F1-score (class 0), F1-score (class 1) and accuracy.  721 \n 722 \nEvaluation Metrics 723 \n 724 \nPrimary reported metrics were weighted -F1, F1 class0, F1 class1 and accuracy with 95% CIs obtained via 725 \nbootstrapping as described above. To characterize the benefit of frozen FM embeddings over the baseline 726 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nTPM we report paired effect sizes for weighted-F1, F1class0, F1class1 and accuracy, computed paired by seed 727 \nwithin fold and summarized with the same bootstrap (paired 95% CIs). 728 \n 729 \nExplainability  730 \n 731 \nSHAP-based gene attribution (TPM models). 732 \n 733 \nWe used Shapley Additive exPlanations (SHAP) [ 26, 27] as implemented in AutoXAI4Omics to quantify 734 \nper-gene contributions to the positive-class F1 score. We used the TreeSHAP explainer for the AdaBoost 735 \nclassifier, fitted on training folds only and applied to the held -out test split. To focus on disease -driving 736 \nsignal (rather than model error or protective evidence), we collected and aggregated SHAP scores from all 737 \nfolds (n=3) and all runs (n=5) of the AdaBoost classifier only for true-positive test samples (ground-truth 738 \nis ‘disease’ and predicted positive threshold >= 0.5). To allow comparison of SHAP scores across samples, 739 \naggregated true-positive SHAP scores were subject to  a positive-only normalization, where positive SHAP 740 \nscores were normalized per sample by scaling each sample’s positive SHAP vector into a distribution 741 \nsummed to  1(only positive SHAP scores were retained to allow focus on features driving disease 742 \nprediction)  Normalized SHAP values were then averaged across folds and seeds to generate a normalized 743 \naverage SHAP value per feature.  744 \n 745 \nGeneformer attribution (classifier-aligned integrated gradients) 746 \n 747 \nFocusing on our best performing FM Geneformer, in this study we used transcriptomic information as input 748 \ninto the FM to derive embeddings (see Methods: Extraction of FM embeddings, Geneformer) . We then 749 \nused these embeddings to train a downstream classifier (AdaBoost) to predict endometriosis status , but 750 \nnotably this classifier was not fine-tuned jointly with Geneformer. Using explainability techniques, our goal 751 \nwas to trace disease predictions back through the FM encoder to identify which input genes were most 752 \nresponsible for the classification decision from our resultant AdaBoost model. To derive this explainability, 753 \nwe developed an approach that we call classifier-aligned integrated gradients (CA-IG) to attribute the FM 754 \ntoken inputs (gene tokens) to the downstream disease score used for classification  using Integrated 755 \nGradients (IG), a well-established attribution method (65). However, the “classifier-aligned” aspect of our 756 \napproach, computing gradients using SHAP values rather than through the entire model, brings novelty and 757 \nutility, dramatically reducing computational cost while maintaining gene -level information about what 758 \ndrives predictions.   759 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n 760 \nWhen a tokeni sed RNA-seq sample is processed by Geneformer , it produces a matrix of dimension 761 \n𝐿 × 𝐷_𝑚 (where L represents length and D_m is the embedding size). The Geneformer encoder at layer ℓ 762 \n(the penultimate layer in all our experiments) transforms this into hidden states of dimension L×D. We 763 \nfocused on real tokens only, excluding special tokens like CLS, EOS, and PAD, and appl ied a masking 764 \nprocedure to obtain the sample embedding. Therefore, the sample embedding, denoted as 𝑧(𝑋), was created 765 \nby mean-pooling the masked hidden states over all real tokens, weighted by their mask values 𝑚 ∈ {0,1}^𝐿. 766 \nThis produced a D-dimensional vector (768 dimensions in this instance) that represented the entire sample 767 \nin the embedding space. 768 \n 769 \nOur downstream classifier AdaBoost operates on this D -dimensional embedding space to produce a 770 \npositive-class score (the disease probability). Since the classifier is not fine-tuned jointly with Geneformer 771 \nwe used a sample-specific surrogate direction derived from the reali sed SHAP attribution vector of the 772 \nclassifier's positive logit with respect to this embedding. This surrogate direction 𝑤(𝑥) exists in the same 773 \nD-dimensional space as the embedding.  We then defined a scalar objective 𝑠(𝑋) as the inner product 774 \nbetween this SHAP direction and the sample embedding 𝑧(𝑋). This scalar objective is locally faithful to 775 \nthe classifier and enable d efficient backpropagation into the token embeddings without requiring 776 \ndifferentiation through the separate classifier. We call this approach classifier-aligned integrated gradients 777 \n(CA-IG) because it aligns the gradient direction with the classifier's behaviour. 778 \n 779 \nComputationally, this approach converts the vector sensitivity into a single directional derivative along the 780 \nSHAP direction 𝑤(𝑥). Using the linearity of gradients, computing the gradient of this scalar objective with 781 \nrespect to the embedding X can be expressed as a weighted sum of per -dimension gradients, where each 782 \ndimension is weighted by the corresponding SHAP weight. This means a single backward pass per IG step 783 \nyields the same weighted gradient that a user  would obtain by computing per -dimension gradients 784 \nseparately and weighting them afterward. This reduces the computational cost from 𝑂(𝑆 · 𝐷) to 𝑂(𝑆) 785 \nbackpropagations for S integration steps across embedding dimension D. To prevent data leakage, the 786 \nSHAP direction 𝑤(𝑥) is computed using only the training fold, then applied to held-out test samples. 787 \n 788 \nIG is a path-based attribution method that accumulates gradients along a straight-line path from a baseline 789 \nto the actual input. We used a PAD baseline B of dimension 𝐿 × 𝐷_𝑚, where each position was set to the 790 \nlearned pad token vector. The IG from baseline B to input X were calculated by integrating the partial 791 \nderivative of the objective with respect to the input along the path from B to X, parameterised by α from 0 792 \nto 1. This continuous integral was approximated using a Riemann sum with S=32 discrete steps along the 793 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \npath. At each step, we evaluated the gradient at a point that interpolated between the baseline and the actual 794 \ninput. Token-level scores were obtained by summing these IG across all channels (the D_m embedding 795 \ndimensions) for each token position, resulting in a score 𝑎_𝑖^ + (𝑥) for each gene token i. These token 796 \nscores were then mapped to genes by removing special tokens and using the tokeniser's gene dictionary to 797 \nform gene-level attributions 𝑔^ ± (𝑥). 798 \n 799 \nAkin to the SHAP analysis, we focused on disease-driving signals by reporting summaries that emphasised 800 \npositive (disease-promoting) contributions. For each sample x, we normali sed the attributions per sample 801 \nbefore aggregating. We created a positive-only summary by taking the maximum of each gene's attribution 802 \nand zero, then normali sing these positive values to sum to 1 (with a small constant ε added to the 803 \ndenominator for numerical stability). This g ave us ĝ_𝑖^ + (𝑥), representing gene i's relative positive 804 \ncontribution to the disease prediction for sample x.  We also compute d a signed summary ĝ_𝑖(𝑥)\tthat 805 \npreserved the direction of effects by normali sing using the sum of absolute values across all genes, again 806 \nwith the stabilising constant ε. We attributed only true-positive patients—specifically, test samples where 807 \nthe ground truth y=1 (disease present) and the predicted probability exceeded 0.5. For outer CV fold k and 808 \nrepeat r, we identif ied the set 𝑇_𝑘, 𝑟 of true -positive test samples. We then average d the normali sed 809 \nattributions within each run, within each fold, and finally across all folds. 810 \n 811 \nThe final importance score 𝐼𝐺𝐼𝑚𝑝_𝑖^ + for gene i was calculated by averaging over K folds, R repeats, and 812 \nall true-positive samples within each fold -repeat combination, using the positive -only normalised scores. 813 \nSimilarly, 𝐼𝐺𝐼𝑚𝑝_𝑖^ ± used the signed normalised scores to capture directional effects.  This multi-level 814 \naveraging produced stable, generali sable gene importance scores that reflect consistent patterns across 815 \ndifferent data splits and random initialisations, minimising the influence of any single model or sample. To 816 \nour knowledge, CA -IG is the first attribution protocol tailored to frozen RNA foundation -model 817 \nembeddings with a separate downstream classifier, and it reduces the per-sample IG complexity from O(S 818 \n·D) to O(S) while remaining CPU-efficient. We report IGImp+ in the main text and provide IGImp± in the 819 \nSupplementary material.  820 \n 821 \nGSEA from attribution ranks 822 \n 823 \nTo identify enriched biological pathways associated with  disease we performed gene set enrichment 824 \nanalysis (GSEA) using pre-ranked gene lists derived from feature importance scores. Four separate analyses 825 \nwere performed with results from the following approaches: a) within-cohort TPM model and associated 826 \ngene SHAP scores; b) cross -cohort TPM model and associated gene SHAP scores; c) within-cohort FM 827 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nembedding model (Geneformer) and associated gene attribution scores from CA -IG and d) cross -cohort 828 \nFM embedding model (Geneformer) and associated gene attribution scores from CA -IG. For TPM-based 829 \nanalyses gene importance scores were obtained from SHAP values of the disease classification model whilst 830 \nFM GSEA analyses were run with gene attribution scores calculated from the classifier-aligned integrated 831 \ngradients method. The gene universe was defined as all genes present in the input TPM matrix. For each of 832 \nthe four analyses, gene importance scores were mapped to the gene universe and genes were ranks in 833 \ndescending order of their signed importance. Pre-ranked GSEA was performed using GSEApy  (66) with 834 \n10,000 permutations, testing using gene se ts from the Enrichr GO Biological Processes  library (67). 835 \nPathways significantly enriched at either end of the ranked gene list were identified at FDR < 0.0 5, with 836 \npositive normalized enrichment score (NES) indicating enrichment among disease -associated genes and 837 \nnegative normalized enrichment score (NES) indicating enrichment among healthy -associated genes. 838 \nParameters used are as follows : min_size= 15, max_size= 500, permutation_num= 10,000, Benjamini –839 \nHochberg FDR control (seed= 42). For each term we report the NES, nominal P value, and FDR-corrected 840 \nP-value. “Leading-edge fraction” denotes the leading-edge size divided by the gene-set size. 841 \n 842 \nKnowledge graph creation and LLM inference 843 \n 844 \nTo identify unexplored relationships between our most explainable genes and endometriosis in a sparsely 845 \ncharacterised research area, we applied a literature -driven knowledge graph and large language model -846 \nbased inference method to systematically integrate and reason over published evidence. This yielded a 847 \nliterature-grounded representation of gene -endometriosis relationships suitable for potential hypothesis 848 \ngeneration. Each candidate gene was treated as an independent query entity and used to retrieve all PubMed-849 \nindexed articles in which it was mentioned, using open -source software DeepSearch 850 \n(https://ds4sd.github.io/). From the resulting corpus, we extracted paragraph-level text segments containing 851 \nco-occurrences of each entity and the endometriosis. These texts were processed using named entity 852 \nrecognition and relation extraction to identify biomedical entities such as genes, diseases, proteins, 853 \npathways, and used to identify the semantic relationships amongst them. All entities and relationships were 854 \nrepresented as entity-entity-relation triples, and integrated into a unified knowledge graph, with literature 855 \nprovenance retained at document and paragraph levels. All inferred associations were constrained by the 856 \nunderlying graph structure and traceable to one or more publications. To enable LLM-based graph-enabled 857 \ninference to identify higher -order or implicit relationships supported by the literature and interpret the 858 \nshared relational patterns in the compact KG, we supplied a compact JSON summary to a locally deployed 859 \nLLM (GPT-20b) executed using Ollama inference engine. We also calculated relation-style similarity and 860 \nshared-entity overlap for all our gene pairs (see below) and provided that to the LLM to aid inference.  861 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n 862 \nRelation-style similarity (weighted Jaccard similarity) 863 \n 864 \nTo quantify similarity of connected entities and relations between the five important genes in the knowledge 865 \ngraph, we calculated a weighted Jaccard similarity for each pair of important genes, quantifying the 866 \nsimilarity of their relation-type profiles across all shared entities. For each pair of important genes and for 867 \neach entity common to both genes, we calculated an intersection weight by summing the minimum 868 \nfrequency observed for each relation type, as follows: 869 \n 870 \n𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑖𝑜𝑛 = J\n!\n𝑚𝑖𝑛K𝑐\"(𝑟), 𝑐#(𝑟)L 871 \n 872 \nNext, for that specific entity we captured the union weight which measures total combined relation usage, 873 \nso taking the maximum of the two counts as follows: 874 \n 875 \n𝑢𝑛𝑖𝑜𝑛 = J\n!\n𝑚𝑎𝑥K𝑐\"(𝑟), 𝑐#(𝑟)L 876 \n 877 \nWe then calculate the Jaccard similarity (intersection/union) with 1 indicating that the two genes share 878 \nexactly the same relation types (for shared entities) with same frequencies, and 0 indicating the two genes 879 \nshare no common relation types for any shared entities.  880 \n 881 \nNormalized shared entity overlap 882 \n 883 \nFor each gene pair, to account for substantial variation in the number of triples per gene due to varied 884 \nrepresentation in the literature, we also calculated shared entity overlap directionally as the proportion of 885 \nshared entities from gene1 that were also present in gene2, and reciprocally from gene2 to gene 1. This 886 \nnormalization allows overlap to be evaluated relative to each genes individual representation in the corpus. 887 \n 888 \nShared-entity overlap (Jaccard similarity) 889 \n 890 \nTo complement the relation -style similarity we also calculate a shared-entity Jaccard to quantify shared 891 \nentities between each gene pair. For each gene pair, this is defined as the proportion of shared (overlapping) 892 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nentities relative to the total number of distinct entities linked to each gene.  For two genes 𝑔\" and 𝑔#, the 893 \nshared-entity jaccard is calculated as follows: 894 \n 895 \n𝐽$%&'&'$( (𝑔\", 𝑔#) =\n|𝐸(𝑔\") ∩ 𝐸(𝑔#)|\n|𝐸(𝑔\") ∪ 𝐸(𝑔#)| 896 \n 897 \n 898 \n 899 \nDeclarations  900 \n 901 \nEthics approval and consent to participate 902 \nNot Applicable – all data used in this study was publicly available from the Gene Expression Omnibus 903 \n(GEO)  904 \n  905 \nConsent for publication 906 \nNot Applicable 907 \nAvailability of data and materials  908 \nAll datasets used in this study are bulk RNA -seq data that is publicly available, Table 1 gives a 909 \nsummary of the cohorts used in this study , detailing Gene Expression Omnibus (GEO) ID’s. The 910 \nfoundation models and AutoXAI4Omics used in this study are all open source  (see citations in 911 \ntext). However, we describe a new interpretability method for the FM embeddings that we call 912 \nCA-IG. We will release this code as an open source repo (at https://github.com/IBM/) upon 913 \nacceptance of publication.  914 \n 915 \nCompeting Interests 916 \nThe authors declare that they have no competing interests. Commercial affiliations do not alter 917 \nadherence of the authors to the journal policies on sharing data and materials. 918 \n 919 \nFunding  920 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \nThis work was supported by the Hartree National Centre for Digital Innovation, a collaboration 921 \nbetween STFC and IBM funded by UKRI. 922 \n 923 \nAuthors’ Contributions 924 \nLJG & JK designed the study, LJG/JK/RT/MA/SC/AE contributed towards workflow conceptualization, 925 \nLJG, JK & NM performed the AI technical analysis, NM, JK & LJG performed the KG technical work with 926 \ncontributions from JB, NC created the CA-IG methodology supervised by LJG & JK, JK, LJG & NM wrote 927 \nand edited the manuscript, MM contributed to BMFM embedding extraction, transcriptomic data collection 928 \nand curation, all authors reviewed code and/or manuscript. 929 \n 930 \n 931 \nAcknowledgements 932 \nThis work was supported by the Hartree National Centre for Digital Innovation, a collaboration 933 \nbetween STFC and IBM. 934 \n 935 \n 936 \n 937 \n 938 \n 939 \nReferences 940 \n 941 \n1. Rogers PA, D’Hooghe TM, Fazleabas A, Gargett CE, Giudice LC, Montgomery GW, et al. 942 \nPriorities for endometriosis research: recommendations from an international consensus 943 \nworkshop. Reproductive sciences. 2009;16(4):335–46.  944 \n2. Cox E. Dismissed, ignored and belittled”: the long road to endometriosis diagnosis in the 945 \nUK. Endometriosis UK. 2024;  946 \n3. Zondervan KT, Becker CM, Missmer SA. Endometriosis. New England Journal of Medicine. 947 \n2020;382(13).  948 \n4. Gete DG, Doust J, Mortlock S, Montgomery G, Mishra GD. Associations between 949 \nendometriosis and common symptoms: findings from the Australian Longitudinal Study on 950 \nWomen’s Health. American Journal of Obstetrics and Gynecology. 2023 Nov;229(5):536.e1-951 \n536.e20.  952 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n5. Meuleman C, Vandenabeele B, Fieuws S, Spiessens C, Timmerman D, D’Hooghe T. High 953 \nprevalence of endometriosis in infertile women with normal ovulation and normospermic 954 \npartners. Fertility and sterility. 2009;92(1):68–74.  955 \n6. Simoens S, Dunselman G, Dirksen C, Hummelshoj L, Bokor A, Brandes I, et al. The burden 956 \nof endometriosis: costs and quality of life of women with endometriosis and treated in 957 \nreferral centres. Human reproduction. 2012;27(5):1292–9.  958 \n7. Nisenblat V, Bossuyt PM, Farquhar C, Johnson N, Hull ML. Imaging modalities for the non-959 \ninvasive diagnosis of endometriosis. Cochrane Database of Systematic Reviews. 2016;(2).  960 \n8. Nisenblat V, Bossuyt PM, Shaikh R, Farquhar C, Jordan V, Scheffers CS, et al. Blood 961 \nbiomarkers for the non-invasive diagnosis of endometriosis. Cochrane Database of 962 \nSystematic Reviews. 1996;2016(5).  963 \n9. Moustafa S, Burn M, Mamillapalli R, Nematian S, Flores V, Taylor HS. Accurate diagnosis 964 \nof endometriosis using serum microRNAs. American journal of obstetrics and gynecology. 965 \n2020;223(4):557-e1.  966 \n10. Zhang H, Zhang H, Yang H, Shuid AN, Sandai D, Chen X. Machine learning-based 967 \nintegrated identification of predictive combined diagnostic biomarkers for endometriosis. 968 \nFront Genet. 2023;14:1290036.  969 \n11. Kucukakcali Z, Akbulut S, Colak C. Prediction of genomic biomarkers for endometriosis 970 \nusing the transcriptomic dataset. World Journal of Clinical Cases. 2025;13(20).  971 \n12. Akter S, Xu D, Nagel SC, Bromfield JJ, Pelch KE, Wilshire GB, et al. GenomeForest: An 972 \nEnsemble Machine Learning Classifier for Endometriosis. AMIA Jt Summits Transl Sci 973 \nProc. 2020;2020:33–42.  974 \n13. Bendifallah S, Dabi Y, Suisse S, Jornea L, Bouteiller D, Touboul C, et al. MicroRNome 975 \nanalysis generates a blood-based signature for endometriosis. Scientific Reports. 976 \n2022;12(1):4051.  977 \n14. Shi S, Huang C, Tang X, Liu H, Feng W, Chen C. Identification and verification of 978 \ndiagnostic biomarkers for deep infiltrating endometriosis based on machine learning 979 \nalgorithms. Journal of Biological Engineering. 2024;18(1):70.  980 \n15. Theodoris CV, Xiao L, Chopra A, Chaffin MD, Al Sayed ZR, Hill MC, et al. Transfer 981 \nlearning enables predictions in network biology. Nature. 2023;618(7965):616–24.  982 \n16. Cui H, Wang C, Maan H, Pang K, Luo F, Duan N, et al. scGPT: toward building a 983 \nfoundation model for single-cell multi-omics using generative AI. Nature methods. 984 \n2024;21(8):1470–80.  985 \n17. Hao M, Gong J, Zeng X, Liu C, Guo Y, Cheng X, et al. Large scale foundation model on 986 \nsingle-cell transcriptomics. bioRxiv. 2023;  987 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n18. Gélard M, Richard G, Pierrot T, Cournède PH. BulkRNABert: Cancer prognosis from bulk 988 \nRNA-seq based language models. bioRxiv. 2024;2024–06.  989 \n19. Dandala B, Danziger MM, Barkan E, Biswas T, Gurev V, Hu J, et al. BMFM-RNA: An 990 \nOpen Framework for Building and Evaluating Transcriptomic Foundation Models. arXiv 991 \npreprint arXiv:250614861. 2025;  992 \n20. Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene expression and 993 \nhybridization array data repository. Nucleic acids research. 2002;30(1):207–10.  994 \n21. Lundberg SM, Lee SI. A unified approach to interpreting model predictions. Advances in 995 \nneural information processing systems. 2017;30.  996 \n22. Djos A, Martinsson T, Kogner P, Carén H. The RASSF gene family members RASSF5, 997 \nRASSF6 and RASSF7 show frequent DNA methylation in neuroblastoma. Mol Cancer. 2012 998 \nDec;11(1):40.  999 \n23. Sangsin A, Kuptanon C, Srichomthong C, Pongpanich M, Suphapeetiporn K, Shotelersuk V. 1000 \nTwo novel compound heterozygous BMP1 mutations in a patient with osteogenesis 1001 \nimperfecta: a case report. BMC Med Genet. 2017 Dec;18(1):25.  1002 \n24. Xu Z, Chen X, Song L, Yuan F, Yan Y. Matrix Remodeling-Associated Protein 8 as a Novel 1003 \nIndicator Contributing to Glioma Immune Response by Regulating Ferroptosis. Front 1004 \nImmunol. 2022 Feb 24;13:834595.  1005 \n25. Zhang M, Cai X, Liu J, Zhou J, Shi Q, Jiang Y, et al. A novel lncRNA lncSAMD11-1: 1 1006 \ninteracts with PIP4K2A to promote endometrial decidualization by stabilizing FoxO1 1007 \nnuclear localization. The International Journal of Biochemistry & Cell Biology. 2022 1008 \nOct;151:106280.  1009 \n26. Greenwood E, Maisel S, Ebertz D, Russ A, Pandey R, Schroeder J. Llgl1 prevents 1010 \nmetaplastic survival driven by epidermal growth factor dependent migration. Oncotarget. 1011 \n2016 Sep 20;7(38):60776–92.  1012 \n27. Al-Hetty HRAK, Jabbar AD, Eremin VF, Jabbar AM, Jalil AT, Al-Dulimi AG, et al. The 1013 \nrole of endoplasmic reticulum stress in endometriosis. Cell Stress and Chaperones. 2023 1014 \nMar;28(2):145–50.  1015 \n28. Huang E, Zhang L, Lou J, Wang X, Chen L. Identified endoplasmic reticulum stress-related 1016 \nmolecular cluster and immune characterization in endometriosis. Sci Rep. 2025 Nov 1017 \n4;15(1):38553.  1018 \n29. Choi J, Jo M, Lee E, Lee DY, Choi D. Dienogest regulates apoptosis, proliferation, and 1019 \ninvasiveness of endometriotic cyst stromal cells via endoplasmic reticulum stress induction. 1020 \nMolecular Human Reproduction. 2020 Jan 1;26(1):30–9.  1021 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n30. Maytin EV, Ubeda M, Lin JC, Habener JF. Stress-Inducible Transcription Factor 1022 \nCHOP/gadd153 Induces Apoptosis in Mammalian Cells via p38 Kinase-Dependent and -1023 \nIndependent Mechanisms. Experimental Cell Research. 2001 Jul;267(2):193–204.  1024 \n31. Ghosh AP, Klocke BJ, Ballestas ME, Roth KA. CHOP Potentially Co-Operates with 1025 \nFOXO3a in Neuronal Cells to Regulate PUMA and BIM Expression in Response to ER 1026 \nStress. El-Rifai W, editor. PLoS ONE. 2012 Jun 28;7(6):e39586.  1027 \n32. Marola OJ, Syc-Mazurek SB, Libby RT. DDIT3 (CHOP) contributes to retinal ganglion cell 1028 \nsomal loss but not axonal degeneration in DBA/2J mice. Cell Death Discov. 2019 Oct 1029 \n10;5(1):140.  1030 \n33. Yu X, Li W, Sun S, Li J. DDIT3 is associated with breast cancer prognosis and immune 1031 \nmicroenvironment: an integrative bioinformatic and immunohistochemical analysis. J 1032 \nCancer. 2024;15(12):3873–89.  1033 \n34. Wang S, Hou P, Pan W, He W, He DC, Wang H, et al. DDIT3 Targets Innate Immunity via 1034 \nthe DDIT3-OTUD1-MAVS Pathway To Promote Bovine Viral Diarrhea Virus Replication. 1035 \nSandri-Goldin RM, editor. J Virol. 2021 Feb 24;95(6):e02351-20.  1036 \n35. Wang S, Hou P, Pan W, He W, He DC, Wang H, et al. DDIT3 Targets Innate Immunity via 1037 \nthe DDIT3-OTUD1-MAVS Pathway To Promote Bovine Viral Diarrhea Virus Replication. 1038 \nSandri-Goldin RM, editor. J Virol. 2021 Feb 24;95(6):e02351-20.  1039 \n36. Söderman J, Berglind L, Almer S. Gene Expression-Genotype Analysis Implicates GSDMA , 1040 \nGSDMB , and LRRC3C as Contributors to Inflammatory Bowel Disease Susceptibility. 1041 \nBioMed Research International. 2015;2015:1–10.  1042 \n37. Hodzic D, Kong C, Wainszelbaum MJ, Charron AJ, Su X, Stahl PD. TBC1D3, a hominoid 1043 \noncoprotein, is encoded by a cluster of paralogues located on chromosome 17q12. 1044 \nGenomics. 2006 Dec;88(6):731–6.  1045 \n38. Wang B, Chen D, Hua H. TBC1D3 family is a prognostic biomarker and correlates with 1046 \nimmune infiltration in kidney renal clear cell carcinoma. Molecular Therapy - Oncolytics. 1047 \n2021 Sep;22:528–38.  1048 \n39. Australian Cancer Study, Australian Ovarian Cancer Study Group, Pharoah PDP, Tsai YY, 1049 \nRamus SJ, Phelan CM, et al. GWAS meta-analysis and replication identifies three new 1050 \nsusceptibility loci for ovarian cancer. Nat Genet. 2013 Apr;45(4):362–70.  1051 \n40. Hodzic D, Kong C, Wainszelbaum MJ, Charron AJ, Su X, Stahl PD. TBC1D3, a hominoid 1052 \noncoprotein, is encoded by a cluster of paralogues located on chromosome 17q12. 1053 \nGenomics. 2006 Dec;88(6):731–6.  1054 \n41. Guitart X, Porubsky D, Yoo D, Dougherty ML, Dishuck PC, Munson KM, et al. Independent 1055 \nexpansion, selection and hypervariability of the TBC1D3 gene family in humans [Internet]. 1056 \nGenetics; 2024 [cited 2025 Dec 15]. Available from: 1057 \nhttp://biorxiv.org/lookup/doi/10.1101/2024.03.12.584650 1058 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n42. Shibel R, Sarfstein R, Nagaraj K, Lapkina-Gendler L, Laron Z, Dixit M, et al. The Olfactory 1059 \nReceptor Gene Product, OR5H2, Modulates Endometrial Cancer Cells Proliferation via 1060 \nInteraction with the IGF1 Signaling Pathway. Cells. 2021 Jun 12;10(6):1483.  1061 \n43. Qin S, Dorschner RA, Masini I, Lavoie‐Gagne O, Stahl PD, Costantini TW, et al. TBC1D3 1062 \nregulates the payload and biological activity of extracellular vesicles that mediate tissue 1063 \nrepair. The FASEB Journal. 2019 May;33(5):6129–39.  1064 \n44. Wagner M, Hicks C, El-Omar E, Combes V, El-Assaad F. The Critical Role of Host and 1065 \nBacterial Extracellular Vesicles in Endometriosis. Biomedicines. 2024 Nov 12;12(11):2585.  1066 \n45. Chen S, Birk DE. The regulatory roles of small leucine‐rich proteoglycans in extracellular 1067 \nmatrix assembly. The FEBS Journal. 2013 May;280(10):2120–37.  1068 \n46. Woo J, Kwon SK, Kim E. The NGL family of leucine-rich repeat-containing synaptic 1069 \nadhesion molecules. Molecular and Cellular Neuroscience. 2009 Aug;42(1):1–10.  1070 \n47. Adachi M, Nasu K, Tsuno A, Yuge A, Kawano Y, Narahara H. Attachment to extracellular 1071 \nmatrices is enhanced in human endometriotic stromal cells: a possible mechanism underlying 1072 \nthe pathogenesis of endometriosis. European Journal of Obstetrics & Gynecology and 1073 \nReproductive Biology. 2011 Mar;155(1):85–8.  1074 \n48. Mariadas H, Chen JH, Chen KH. The Molecular and Cellular Mechanisms of Endometriosis: 1075 \nFrom Basic Pathophysiology to Clinical Implications. IJMS. 2025 Mar 10;26(6):2458.  1076 \n49. Miller JE, Ahn SH, Marks RM, Monsanto SP, Fazleabas AT, Koti M, et al. IL-17A 1077 \nModulates Peritoneal Macrophage Recruitment and M2 Polarization in Endometriosis. Front 1078 \nImmunol. 2020 Feb 14;11:108.  1079 \n50. Shi JL, Zheng ZM, Chen M, Shen HH, Li MQ, Shao J. IL-17: an important pathogenic factor 1080 \nin endometriosis. Int J Med Sci. 2022;19(4):769–78.  1081 \n51. Ahn SH, Edwards AK, Singh SS, Young SL, Lessey BA, Tayade C. IL-17A Contributes to 1082 \nthe Pathogenesis of Endometriosis by Triggering Proinflammatory Cytokines and 1083 \nAngiogenic Growth Factors. The Journal of Immunology. 2015 Sep 15;195(6):2591–600.  1084 \n52. Zhang X, Xu H, Lin J, Qian Y, Deng L. Peritoneal fluid concentrations of interleukin‐17 1085 \ncorrelate with the severity of endometriosis and infertility of this disorder. BJOG. 2005 1086 \nAug;112(8):1153–5.  1087 \n53. Delbandi AA, Mahmoudi M, Shervin A, Farhangnia P, Mohammadi T, Zarnani AH. 1088 \nIncreased circulating T helper 17 (TH17) cells and endometrial tissue IL-17-producing cells 1089 \nin patients with endometriosis compared with non-endometriotic subjects. Reproductive 1090 \nBiology. 2025 Jun;25(2):101019.  1091 \n54. Hirata T, Osuga Y, Takamura M, Saito A, Hasegawa A, Koga K, et al. Interleukin-17F 1092 \nincreases the secretion of interleukin-8 and the expression of cyclooxygenase 2 in 1093 \nendometriosis. Fertility and Sterility. 2011 Jul;96(1):113–7.  1094 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint \n\n   \n \n   \n \n55. Podgaec S, Abrao MS, Dias JA, Rizzo LV, De Oliveira RM, Baracat EC. Endometriosis: an 1095 \ninflammatory disease with a Th2 immune response component. Human Reproduction. 2007 1096 \nMay 1;22(5):1373–9.  1097 \n56. Wu MH, Yang BC, Hsu CC, Lee YC, Huang KE. The expression of soluble intercellular 1098 \nadhesion molecule-1 in endometriosis. Fertility and Sterility. 1998 Dec;70(6):1139–42.  1099 \n57. Keenan JA, Chen TT, Chadwell NL, Torry DS, Caudle MR. Interferon‐Gamma (IFN‐γ) and 1100 \nInterleukin‐6 (IL‐6) in Peritoneal Fluid and Macrophage‐Conditioned Media of Women With 1101 \nEndometriosis. American J Rep Immunol. 1994 Oct;32(3):180–3.  1102 \n58. Al-Hetty HRAK, Jabbar AD, Eremin VF, Jabbar AM, Jalil AT, Al-Dulimi AG, et al. The 1103 \nrole of endoplasmic reticulum stress in endometriosis. Cell Stress and Chaperones. 2023 1104 \nMar;28(2):145–50.  1105 \n59. Ekiz-Yilmaz T, Isildar B, Gezer A, Kankaya D, Cansiz-Ersoz C, Kayisli UA, et al. The role 1106 \nof unfolded protein response in the pathogenesis of endometriosis: contribution of peritoneal 1107 \nfluid. Reproductive BioMedicine Online. 2021 Jan;42(1):1–15.  1108 \n60. Tellez-Quijorna C, Juan-Lopez A, Eritja N, Llobet-Navas D, Devis-Jauregui L. Dysregulated 1109 \nautophagy in endometriosis: molecular mechanisms, controversies, and clinical implications. 1110 \nReprod Biol Endocrinol. 2025 Dec 13;24(1):8.  1111 \n61. Mariadas H, Chen JH, Chen KH. The Molecular and Cellular Mechanisms of Endometriosis: 1112 \nFrom Basic Pathophysiology to Clinical Implications. IJMS. 2025 Mar 10;26(6):2458.  1113 \n62. Hao M, Gong J, Zeng X, Liu C, Guo Y, Cheng X, et al. Large scale foundation model on 1114 \nsingle-cell transcriptomics. bioRxiv. 2023;  1115 \n63. Strudwick J, Gardiner LJ, Denning-James K, Haiminen N, Evans A, Kelly J, et al. 1116 \nAutoXAI4Omics: an automated explainable AI tool for omics and tabular data. bioRxiv 1117 \n[Internet]. 2024; Available from: 1118 \nhttps://www.biorxiv.org/content/early/2024/03/28/2024.03.25.586460 1119 \n64. Freund Y, Schapire RE. A decision-theoretic generalization of on-line learning and an 1120 \napplication to boosting. Journal of Computer and System Sciences. 1997;55(1):119–39.  1121 \n65. Sundararajan M, Taly A, Yan Q. Axiomatic attribution for deep networks. In: International 1122 \nconference on machine learning. PMLR; 2017. p. 3319–28.  1123 \n66. Fang Z, Liu X, Peltz G. GSEApy: a comprehensive package for performing gene set 1124 \nenrichment analysis in Python. Lu Z, editor. Bioinformatics. 2023 Jan 1;39(1):btac757.  1125 \n67. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a 1126 \ncomprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 1127 \n2016 Jul 8;44(W1):W90–7.  1128 \n 1129 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint","source_license":"CC0","license_restricted":false}