RNA foundation models enable generalizable endometriosis disease classification and stable gene-level interpretation

16 17 Endometriosis is a chronic inflammatory condition with significant diagnostic delays impacting one in ten 18 reproductive age women worldwide. While machine learning (ML) models trained on transcriptomic data 19 show promise for disease prediction, limited generalizability across independent patient cohorts has 20 hindered clinical translation. Foundations models (FMs) pretrained on large-scale transcriptomic data offer 21 promise to learn transferrable, biologically meaningful representations that could support cross -cohort 22 predictions. We assembled a 12 -cohort bulk RNA -seq benchmark (334 samples) and developed a 23 computationally efficient pipeline to test whether FMs improve endometriosis classification, an approach 24 not previously applied to this disease. Using AutoXAI4Omics with cohort-aware validation, we compared 25 embeddings derived from five state-of-the-art RNA FMs against TPM baselines. In cross-cohort prediction, 26 FM embeddings significantly improved performance, achieving a weighted F1 -score of 0.83 vs. 0.68 for 27 the baseline. To allow gene-level interpretation of FM embedding models, we introduce classified-aligned 28 integrated gradients (CA -IG), an interpretability approach aligning gene -level attributions to the 29 downstream classifier without end-to-end finetuning. CA-IG revealed a conserved set of predictive genes 30 from FM embeddings across cohort-validation regimes, contrasting with unstable baseline explainability, 31 suggesting that FM embeddings prioritized transferable disease-related signal over cohort-specific effects. 32 These genes include novel candidates that converge on biologically plausible pathways for endometriosis. 33 34 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint

35 Foundation Models, Machine learning, RNA, Transcriptomics, Endometriosis, disease classification, 36 explainability 37 38

39 40 Endometriosis remains a diagnostic and therapeutic challenge, affecting one in ten women of reproductive 41 age and imposing a significant personal and economic burden. Definitive d iagnosis still hinges on 42 laparoscopy and biopsy with expected diagnosis times of approximately 9 years in the UK (1–3). The 43 disease is characterized by ectopic endometrial -like tissue and typically presents with chronic pelvic or 44 abdominal pain , dysmenorrhoea, dyspareunia, fatigue, and gastrointestinal or urinary symptoms (4). 45 Endometriosis is also a leading cause of subfertility/infertility, with a prevalence of 30-50% among women 46 seeking infertility care (5). Health-related quality of life is markedly reduced, with increased anxiety and 47 depression; indirect and direct costs are considerable, and in the UK alone the annual economic burden has 48 been estimated at £8.2 billion (6). Clinical assessment is further complicated by phenotypic heterogeneity 49 and a poor correlation between lesion burden and symptom severity, while imaging is often inconclusive 50 for superficial disease and reliable blood or tissue biomarkers remain lacking (7–9). Hindering diagnostic 51 development is the fact that precise molecular drivers of endometriosis physiology are not well defined or 52 understood. As such, there is an increasing focus and need to interrogate high -throughput omics data like 53 transcriptomics, which offer the ability to study genome-wide expression patterns and reveal dysregulated 54 pathways and genes across different patient groups, lesion types and tissue/fluid types. Insights from these 55 types of studies offer potential for improved diagnostics, biomarker identification, disease subtype 56 classification, increased mechanistic understanding of pathogenesis and personalized medicine approaches. 57 This has motivated the application of AI to predict or model endometriosis from high-throughput expression 58 data, including microarrays and RNA -seq, for various applications such as disease classification and 59 biomarker discovery (9–13). 60 61 Akter et al. used an ensemble-based approach for endometriosis prediction on single-centre bulk RNA-seq 62 data (n=38) and DNA-methylation data, and reported near-perfect internal performance (F1 ~ 0.97) using 63 leave-one-out cross- validation (CV) (12). Revisiting the same cohort, another study used elastic-net feature 64 selection with boosted and bagged trees and achieved performance around F1~0.86 (5 -fold CV) (11). On 65 microarrays, Zhang et al. (2023) proposed a five-gene panel (FOS, EPHX1, DLGAP5, PCSK5, ADAT1) for 66 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint endometriosis classification and achieved a AUC of 0.836 on a held-out test set and > 0.78 on two external 67 validations (10). A parallel line of work has pursued circulating microRNAs as non-invasive liquid-biopsy 68 markers of endometriosis, showing a six -miRNA panel exhibited a high accuracy for endometriosis 69 prediction (AUC=0.94) in an independent 48 -patient validation cohort (9). More recently, Bendifallah et 70 al. profiled the circulating miRNome in a prospective 200 -patient study and derived an 86 -miRNA 71 signature with excellent internal discrimination (AUC~0.98 under cross -validation), but without multi -72 centre external replication (13). As the authors note, the broader blood -biomarker literature has been 73 constrained by small cohorts and methodological heterogeneity, yielding few stable, clinically adopted 74 markers. This finding was enforced by Shi et al. who demonstrated that subtype-specific tissue markers for 75 deep infiltrating endometriosis were sensitive to cohort composition (14). 76 77 While these studies revealed critical biological insights and significantly advanced endometriosis disease 78 understanding, their reliance on single cohorts has limited the generalizability of their models and 79 associated findings to unseen datasets or cohorts. Moreover, this issue limits the clinical applicability and 80 value of their findings in real -world diagnostic settings. As such, there is a need to develop methods or 81 approaches which can derive more generalized insights for endometriosis that will be applicable across a 82 wider range of cohorts and can underpin development of robust models and derivative disease biomarkers. 83 A potential solution to this limitation may be the use of foundation models (FMs), large neural encoders 84 pretrained on vast unlabeled corpora via self -supervision to learn general -purpose representations that 85 transfer across tasks. Recent work has demonstrated the potential of FMs to learn representations that are 86 highly transferrable across different domains, datasets and tasks. In transcriptomics, RNA FMs such as 87 Geneformer (15), scGPT (16), scFoundation (17), and BulkRNABERT (18) are trained on millions of bulk 88 or single-cell gene expression profiles and have shown utility for downstream problems including cell-type 89 annotation, perturbation prediction, and disease classification . Relative to models trained from scratch, 90 these encoders typically yield increased downstream task performance and greater robustness under 91 distribution shift by providing richer, structure -aware features. This study represents the first systematic 92 evaluation of RNA foundation models for endometriosis prediction. 93 94 In this work we aimed to investigate if use of RNA FMs meaningfully improved generalizability of the 95 resulting models, allowing accurate cross-cohort disease prediction. To do this we developed a sophisticated 96 cross-dataset evaluation pipeline to test the domain generalizability of FM embeddings using five state-of-97 the-art RNA FMs; Geneformer (15), scFoundation (17), scGPT (16), BulkRNABERT (18) and BMFM-98 RNA (19) and a frozen-encoder paradigm (pretrained RNA encoders are kept fixed and used to compute 99 embeddings which were then used as input features for a classical ML task). The resultant pipeline is 100 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint detailed in Figure 1 and consisted of 1) identification of Gene Expression Omnibus (GEO) ids for relevant 101 bulk RNA -seq endometriosis/control datasets and automated download, extraction, amalgamation and 102 normalization of transcriptomic datasets given these pre-defined ids, and LLM-powered metadata curation 103 with human validation to ensure consistent and high -quality target label generation; 2) extraction of 104 embeddings from pre-trained frozen RNA FMs via one forward pass through each model; 3) generation of 105 robust cohort-aware multi-fold strategy to evaluate within vs. cross cohort learning capabilities of using 106 FM embeddings as features compared to baseline normalized (TPM) gene expression counts; 4) an 107 explainable ML workflow to train and test a decision tree classifier to predict endometriosis status (see 108 Methods) and to evaluate performance of our feature sets; 5) classifier -aligned integrated gradients (IG) 109

to calculate gene attribution scores and 6) an LLM -powered knowledge graph (KG) approach to 110 interpret biological significance of the most predictive features. 111 112 113 114 115 Figure 1. End to end pipeline. Pipeline consists of 1) identification of GEO ids for relevant 116 endometriosis/control datasets and LLM -powered metadata curation with human validation to ensure 117 consistent and high-quality target label generation ; 2) extraction of embeddings from pre -trained frozen 118 RNA FMs; 3) generation of robust cohort-aware multi-fold strategy to test within vs. cross cohort learning 119 capabilities of feature sets; 4) an explainable ML workflow to evaluate endometriosis classification 120 performance of our feature sets; 5) classifier-aligned integrated gradients (CA-IG) method to calculate gene 121 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint attribution scores and 6) an LLM -powered knowledge graph (KG) approach to interpret biological 122 significance of the most predictive features. 123 124 Our dataset incorporated a multi-cohort bulk RNA-seq benchmark spanning 12 independent GEO studies 125 (334 samples) , one of the largest collection s of independent cohorts for bulk RNA -seq endometriosis 126 classification to date, and the first multi -cohort evaluation of frozen RNA foundation -model embeddings 127 for this task. This approach allowed our data to be enhanced with the prior biological knowledge of the FM, 128 without the computational expense of pre -training or fine -tuning. For interpretation of the embeddings 129 derived from the FMs, we applied a hybrid classifier-aligned integrated gradients (CA-IG) method, which 130 can produce stable, biologically plausible gene -level attributions across evaluation regimes. To identify 131 unexplored relationships between our most explainable genes and endometriosis in a sparsely characterized 132 research area, we applied a literature -driven knowledge graph and LLM-based inference method to 133 systematically integrate and reason over published evidence, yielding a literature-grounded representation 134 of gene-endometriosis relationships suitable for hypothesis generation. 135 136 This work enabled a robust assessment of the capability of FM -based representations to generalize across 137 heterogenous sources, a persistent challenge situation in computational biology, with often limited 138 availability of large, harmonized datasets. Many studies rely on small sample sizes collected under specific 139 conditions, which commonly leads to individual models with good performance but poor generalizability. 140 This lack of generalizability is of particular concern with understudied complex diseases like endometriosis. 141 By using FM-derived embeddings in this way we aim to address this fundamental issue and understand if 142 representations learnt from pre -trained foundation models capture biologically genuine patterns which 143 transcend dataset-specific biases. 144 145 146

& Discussion 147 148 Novel pipeline design enables evaluation of cross-cohort generalizability 149 150 In this work we developed a comprehensive end-to-end pipeline enabling inference-only use of pre-trained 151 RNA foundation models for reproducible endometriosis disease prediction and gene-level interpretation 152 from bulk RNA-seq data (Fig 1). This pipeline integrates: (i ) automated download and harmonization of 153 GEO RNA-seq studies with LLM -assisted metadata curation and expert validation, (ii) extraction of FM 154 embeddings via a single forward pass through the model, (iii) cohort-aware evaluation using within-cohort 155 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint and cross-cohort strategies, (iv) consistent classical ML training and evaluation using AutoXAI4Omics 156 tookit, (v) explainability using SHAP for TPM models and classifier-aligned integrated gradients (CA-IG) 157 for embedding models and (vi) literature -grounded knowledge-graph generation and interpretation, and 158 pathway analysis. We assembled a multi -cohort benchmark spanning 12 independent Gene Expression 159 Omnibus GEO studies (334 samples ; 17,518 unified genes ) (20), by searching for human studies of 160 endometrial tissue related to endometriosis and restricting results to expression profiling by high -161 throughput sequencing, comprising 259 endometriosis cases and 75 control s (Table 1). This multi -study 162 setting captures clinically realistic heterogeneity across laboratories and protocols and enables an accurate 163 assessment of generalization under cohort shift. 164 165 From these datasets we derived two feature representations for each sample; (i) gene-level log-transformed 166 TPM values, used as a baseline, and (ii) embeddings extracted from state-of-the-art ‘frozen’ RNA FMs. To 167 assess if FM-derived features provide an advantage over baseline features, we evaluated two classification 168 approaches. In approach 1, we used samples from all 12 cohorts to build the training and test sets, simulating 169 a ‘within-cohort’ approach (see Methods). Here we expect elevated performance since we do not ask the 170 model to generalize to a completely unseen cohort, but in doing so we risk models learning patterns distinct 171 to each dataset . For approach 2, we designed test and train sets which consisted of distinct cohorts, 172 simulating a ‘cross-cohort’ approach (see Methods). In this approach we trained the model on samples 173 from specified cohorts and tested exclusively on samples from independent cohorts which were not seen 174 during training. Here we expect a drop -off in performance as the model must generalize across many 175 cohorts. The latter scenario reflects a common challenge in clinical and biological research, where data are 176 typically distributed across many small, independent cohort. By comparing performance across these two 177 strategies, we aimed to determine whether FM-derived embeddings improve robustness and generalization 178 in the presence of cohort heterogeneity, relative to classical gene-level features as a baseline. 179 180 181 Table 1: Summary of the cohorts used in this study. Detailing GEO ID’s, sample counts and positive 182 proportions of patients or samples within each cohort with a positive diagnosis of endometriosis. 183 184 Dataset ID Dataset name (GEO name) Sample count(n) Proportion of disease cases 1 GSE134056 38 0.42 2 GSE202571 17 1 3 GSE205494 12 0.5 4 GSE135485 58 0.93 5 GSE212787 17 0.65 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 6 GSE135640 36 0.5 7 GSE226870--GSE226872 53 1 8 GSE157153 66 1 9 GSE230956 8 0.5 10 GSE171653 11 1 11 GSE190580 7 0.43 12 GSE243157--GSE243158 11 0 185 Within-cohort disease prediction 186 187 We first evaluated the within-cohort performance (where training and test sets include samples from the 188 same cohort, though not the same patient ) for endometriosis disease classification (Table 2). To obtain a 189 robust analysis of model performance variability and to allow statistical comparison, we evaluated 190 performance by performing 3-cross validation folds with 5 random-seed runs per fold per model (15 runs 191 in total per model), and performed hierarchical bootstrapping (10,000 iterations) with re-sampling, allowing 192 calculation (for each metric) of mean, mean effect size , 95% confidence interval (CI) and two-tailed p-193 values (see Methods). 194 195 ML classification of endometriosis disease using the log transformed TPM data (baseline) resulted in a 196 mean 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 accuracy of 0.86 ( Table 2, Fig 2A), indicating good performance for endometriosis classification from 198 multi-cohort transcriptomic data, with more accurate prediction of disease than healthy samples (Table 2). 199 200 201 Table 2. Within-cohort results . Predictive performance metrics for endometriosis disease prediction 202 using our within-cohort strategy. Four model performance metrics are provided (weighted -F1, class0-F1, 203 class1-F1 and accuracy) from classification models built using the baseline TPM features and FM -204 embedding features from five RNA FMs. For GeneFormer the V2 -104M-HC model was used and for 205 BMFM-RNA the 110m -mlm-rda-v2 model was used. For each metric, table reports the estimated mean 206 performance over 10,000 bootstrap replicates and corresponding low (2.5 th) and high (97.5th) values from 207 calculating the 95% confidence interval (CI) obtained via nested bootstrapping which re-sampled folds and 208 seed runs with replacement. All metrics are computed on the held -out test set, and ML models were built 209 and evaluated using AutoXAI4Omics ML pipeline (see Methods), using AdaBoost. 210 211 Model Name Weighted F1 Class 0 – F1 Class 1 – F1 Accuracy .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint (mean [95%CI]) (mean [95%CI]) (mean [95%CI]) (mean [95%CI]) TPM 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] BulkRNABert 0.89 [0.88,0.91] 0.77 [0.75, 0.79] 0.93 [0.92,0.94] 0.89 [ 0.88, 0.91] scFoundation 0.85 [0.81, 0.87] 0.63 [0.55,0.69] 0.91 [0.89,0.92] 0.85 [0.82, 0.88] scGPT 0.85 [0.81, 0.87] 0.65 [0.57,0.72] 0.91 [0.89,0.93] 0.86 [0.82, 0.89] GeneFormer 0.90 [0.86, 0.91] 0.76 [0.69,0.81] 0.93 [0.91, 0.94] 0.89 [0.86, 0.91] BMFM-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] 212 213 Next, we computed the statistical significance of pairwise model performance by calculating two-tailed p-214 values from the bootstrapped samples. Among the five frozen FMs, embeddings extracted from 215 BulkRNABERT and Geneformer surpassed the predictive performance of the TPM baseline classifier in 216 all four metrics with statistical significance (P < 0.05, Supplementary Table S1A), with the weighted F1 217 increasing from 0.86 for the TPM classifier to 0.89 and 0.90 for BulkRNABERT and Geneformer classifiers 218 respectively (Table 2). Prediction of the healthy class (class 0) was also significantly improved using 219 BulkRNABERT and Geneformer embedding classifiers, increasing from a baseline-F1 of 0.69 to F1 scores 220 of 0.77 and 0.76 for BulkRNABERT and Geneformer classifiers respectively (Table 2, Supplementary 221 Table S2A, Fig 2). For the remaining three FM embedding classifiers (scFoundation, scGPT and BMFM -222 RNA) we saw reduced or equal model performance compared to the baseline TPM classifier for all metrics, 223 and for prediction of the healthy class (class 0), the scFoundation embedding classifier performed 224 significantly worse than the TPM model (P < 0.05, Supplementary Tables S2A, S3A and S4A). 225 Surprisingly, our top performing FM-embedding classifier was using embeddings extracted from the single-226 cell pretrained FM Geneformer, whereas we were using bulk RNA-seq data. This suggests that single-cell 227 priors transfer effectively to bulk RNA-seq. 228 229 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 230 231 Figure 2: Within-cohort results. (a) Mean predictive performance metrics for endometriosis using our 232 within-cohort strategy for the baseline TPM features versus using frozen embeddings as from the five RNA 233 FMs as features. Bars show the mean across outer folds and seeds; whiskers denote 95% CIs. Metrics 234 reported: weighted-F1 (primary), F1 for the positive class (class 1; endometriosis), F1 for the negative class 235 (class 0; controls), and accuracy. (b) Effect size relative to the baseline: mean difference in score 236 (embedding model minus baseline) with 95% Cis. All models use identical training and evaluation within 237 the series-aware within-cohort splits 238 239 For each FM-derived feature set classifier, mean effect-size (+ 95% CI) was also calculated which indicated 240 the average difference in performance metrics from the raw baseline transcriptomic feature set model, 241 which we treated as our baseline (Fig. 2B, Supplementary Table S5). Although we observed an increase in 242 F1 score and positive F1 effect sizes for two FM -derived embedding classifiers, increases were moderate 243 in scale with very small effect sizes, and we observed a drop in model performance for the remaining three 244 FM embedding classifiers. Overall, within -cohort evaluation yielded high performance for both feature 245 representation type (TPM and embeddings), with only modest and model -dependent benefits from 246 embeddings (Fig 2B). This pattern is consistent with a setting where cohort -specific structure is shared 247 between test and train data, and in these conditions a well-tuned classical pipeline using gene-level features 248 already performs strongly, with embeddings adding limited advantage. 249 250 FM-derived embeddings substantially improve cross-cohort disease prediction 251 252 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint We next evaluated cross-cohort generalization, where entire cohorts are held out for testing. For this 253 approach, we constructed balanced grouped cohort splits with K=3 outer folds to match the within-cohort 254 approach (see Methods). Individual cohorts were assigned to test folds by a constrained random search 255 (1,000 trials; fixed seed) that (i) enforced a minimum test size (at least 60 samples per fold) and (ii) 256 minimized the across-fold standard deviation of the positive fraction (disease positive samples) and of fold 257 sizes. We report the resultant test/train dataset assignments in Methods. Mean model performance metrics 258 and effect sizes were calculated as previously described (see Methods). 259 260 When using our cross-cohort approach for disease prediction, the baseline classifier (log-transformed TPM) 261 achieved a mean weighted F1-score of 0.68, a marked decline from the F1-score of 0.86 obtained from the 262 within-cohort approach (Table 3). This loss in performance was expected, as the trained model was less 263 able to generalize to samples from unseen cohorts, likely due to differences in sample collection, patient 264 clinical diversity, feature distribution and class-imbalances. All embedding-based classifiers also exhibited 265 a drop in predictive performance when trained and tested using the cross-cohort vs. within-cohort approach, 266 however the degree of reduction was substantially lower in 4/5 of the embedding classifiers compared to 267 the baseline classifier, with drops in mean weighted F1-score of 0.07, 0.07, 0.10 and 0.11 for GeneFormer, 268 scGPT, BMFM-RNA and BulkRNABert models respectively, compared to a drop of 0.18 for the baseline 269 (Table 3). 270 271 272 Table 3. Predictive performance metrics from endometriosis disease prediction (cross-cohort). Four model 273 performance metrics are provided (weighted-F1, class0-F1, class1-F1 and accuracy) from models using the 274 baseline TPM features and FM -embedding features from 5 SOTA transcriptomic FMs. For GeneFormer 275 the V2-104M-HC model was used and for BMFM-RNA the 110m-mlm-rda-v2 model was used. For each 276 metric, table reports the estimated mean performance over 10,000 bootstrap replicates and corresponding 277 low (2.5th) and high (97.5th) values from calculating the 95% confidence interval (CI) obtained via nested 278 bootstrapping which re -sampled folds and seed runs with replacement. All metrics are computed on the 279 held-out test set, and ML models were built and evaluated using AutoXAI4Omics ML pipeline (see 280 Methods), using AdaBoost. 281 282 Model Name Weighted F1 (mean[95%CI]) Class 0 – F1 (mean [95%CI]) Class 1 – F1 (mean [95%CI]) Accuracy (mean [95%CI]) TPM 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] BulkRNABert 0.78 [0.68, 0.90] 0.53 [0.21, 0.73] 0.86 [0.82,0.91] 0.79 [ 0.71, 0.86] .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint scFoundation 0.58 [0.52, 0.63] 0.31 [0.08,0.41] 0.66 [0.55,0.79] 0.57 [0.50,0.66] scGPT 0.78 [0.74, 0.81] 0.50 [0.32,0.61] 0.86 [0.85,0.87] 0.79 [0.77, 0.80] GeneFormer 0.83 [0.77, 0.87] 0.62 [0.44,0.74] 0.89 [0.86, 0.91] 0.83 [0.78,0.87] BMFM-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] 283 284 In contrast to our within -cohort approach, now when using a cross -cohort approach, all FM -embedding 285 classifiers significantly outperformed the TPM -baseline model for disease classification, except for 286 scFoundation ( P < 0.05, Supplementary Table S1B). Geneformer and BMFM -RNA were the top two 287 performing embedding models for all metrics (mean weighted F1-score, class 0 F1-score, class 1 F1-score 288 and mean accuracy) with mean weighted F1-scores of 0.83 and 0.80 respectively (Table 3, Fig 3, P < 0.05). 289 Significant increases in performance for all embedding models apart from scFoundation compared to 290 baseline were also apparent for class 1 F1 -score and accuracy metrics ( P < 0.05), while GeneFormer, 291 BMFM-RNA and scGPT also achieved significantly higher F1-scores for class 0 compared to baseline (P 292 < 0.05, Table 3, Supplementary Table S1B). For the outperforming FMs, we saw an average positive effect 293 size (across all four models) compared to the baseline of 0.11 for weighted F1, class 0 F1-score and class1 294 F1-score and of 0.13 for accuracy (Supplementary Table S6). We saw the greatest performance 295 improvement from Geneformer and BMFM -RNA embeddings, where class 0 F1 increased (compared to 296 TPM baseline) by an average of 0.2 and 0.14 respectively, and the weighted F1 -score increased by an 297 average of 0.15 and 0.12 respectively, reflecting a substantial and significant improvement in model 298 performance (Table 3, Fig 3, Supplementary Table S6). 299 300 301 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 302 303 Figure 3: Cross-cohort results (balanced grouped -cohort validation). (a) Mean performance for the 304 (TPM-Baseline) versus the same pipeline using frozen embeddings (BulkRNABERT, scFoundation, 305 scGPT, Geneformer, BMFM -MLM-RDA). Bars show the mean across outer folds and seeds; whiskers 306 denote 95% CIs. Metrics: weighted -F1 (primary), F1 for the positive class (class 1), F1 for the negative 307 class (class 0), and accuracy. (b) Effect size relative to the baseline: mean score difference (embedding 308 model minus baseline) with 95% CIs; the vertical dashed line at 0 marks parity (values > 0 favour 309 embeddings). In this setting, train and test cohorts are disjoint: each outer fold holds out a group of GEO 310 Series (i.e., of cohorts) for testing, and all model selection is confined to the training side 311 312 Gene-level attribution enables interpretation of baseline TPM model 313 314 We investigated explainability of the within-cohort and cross-cohort baseline TPM classifiers using SHAP 315 (Shapley Additive Explanations) a model -agnostic attribution framework which compute s feature 316 importance scores using game theory (21). This approach allows us to identify the genes most important 317 for disease prediction. Using results from all folds (n=3) and all runs (n=5) per model, for each true positive 318 sample (endometriosis disease positive, correctly predicted as disease) we first extracted SHAP values for 319 all genes and normalized across folds . To focus on disease -driving signal o ur normalization approach 320 included first setting all negative SHAP scores to zero and normalizing the remaining positive samples, 321 resulting in proportionate importance scores for disease-driving genes. Normalized SHAP values were then 322 averaged across models, runs and folds, resulting in one average SHAP value per gene per model. For the 323 within-cohort baseline TPM model (AdaBoost) across all runs and folds there were an average of 75 true 324 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint positive samples from which the SHAP values were normalized and averaged. A total of 51 genes exhibited 325 non-zero SHAP values, with the top 20 disease-associated genes shown in Table 4. For the cross-cohort 326 baseline TPM model (AdaBoost) across all runs and folds there were an average of 111 true positive 327 samples from which the SHAP values were normalized and averaged with a total of 64 genes that exhibited 328 non-zero SHAP values (Table 4). 329 330 331 Table 4. Baseline TPM explainability. Top 20 positive (disease-driving) genes identified via SHAP for 332 TPM baseline model for the within -cohort approach and compared to the cross -cohort approach. 333 Normalized using a positive -only approach, then averaged across runs and folds (HGNC symbols) (Full 334 gene list available in Supplementary Tables S7 (within-cohort) and S8 (cross-cohort)). 335 336 WITHIN COHORT CROSS COHORT GeneID score Gene name GeneID score Gene name EFEMP2 0.17941397 Fibulin 4-encoding gene UCN2 0.0951607 urocortin 2 RASSF6* 0.08950262 Ras Association Domain Family Member 6-encoding gene PCOLCE 0.09124098 procollagen C- endopeptidase enhancer SLC39A13 0.0727763 Solute Carrier Family 39 Member 13 (zinc transporter) BMP1* 0.07026173 bone morphogenetic protein 1 TP53I13 0.064122 Tumor Protein P53 Inducible Protein 13 LLGL1* 0.05474552 LLGL scribble cell polarity complex component 1 BMP1* 0.05455228 Bone Morphogenetic Protein 1- encoding gene MAL 0.05190146 mal, T cell differentiation protein (MAL blood group) LCN2 0.04278986 Lipocalin 2 -encoding gene CYP8B1 0.05036643 cytochrome P450 family 8 subfamily B member 1 STEAP4 0.03793841 STEAP4 Metalloreductase RASSF6* 0.04911975 Ras association domain family member 6 VPS16 0.03228817 Vacuolar Protein Sorting Protein 16-encoding gene TSPAN4 0.04771203 tetraspanin 4 CDC37 0.0305804 Cell Division Cycle 37, HSP90 Cochaperone AGRN 0.03964136 agrin SAMD11* 0.03008721 Cell division cycle 37 homolog SAMD11* 0.03699761 Cell division cycle 37 homolog .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint MXRA8* 0.02810625 Matrix remodeling-associated protein 8 ARFGAP1 0.03572216 ARF GTPase activating protein 1 OAF 0.02796728 Out at first homolog IDO1 0.03315175 indoleamine 2,3- dioxygenase 1 PDK4 0.02608833 pyruvate dehydrogenase kinase 4 CCDC141 0.03049321 coiled-coil domain containing 141 SURF4 0.02398875 surfeit 4 CD276 0.02285808 CD276 molecule HPS1 0.02067466 HPS1 biogenesis of lysosomal organelles complex 3 subunit 1 MLXIPL 0.02278676 MLX interacting protein like TMEM273 0.01807645 transmembrane protein 273 MXRA8* 0.02076158 matrix remodeling associated 8 TAX1BP3 0.01622053 Tax1 binding protein 3 ITGAM 0.01675609 integrin subunit alpha M LLGL1* 0.01517964 LLGL scribble cell polarity complex component 1 CHST14 0.01642441 carbohydrate sulfotransferase 14 ADAMTS2 0.013455 ADAM metallopeptidase with thrombospondin type 1 motif 2 TRAF7 0.01502836 TNF receptor associated factor 7 LYZ 0.01336511 lysozyme TMEM129 0.01500533 transmembrane protein 129, E3 ubiquitin ligase 337 338 We observed little overlap between the TPM-baseline explainability results for the within-cohort vs cross-339 cohort approaches, with only 5/20 of the top 20 explainable genes in the within-cohort approach also within 340 the top 20 for the cross-cohort approach (Table 4). The top four most important genes for the within-cohort 341 and cross-cohort approaches were completely distinct, being derived from different biological pathways 342 (WikiPathways), indicating that different biological signals were being captured by the within and cross-343 cohort approaches (Supplementary Figure 1, Table 4). This variability in the predictive features between 344 models is concerning for researchers using interpretation of ML models to guide biomarker or target 345 discovery for disease and could partially explain our drop in performance for cross-cohort analysis if the 346 predictive features from the within-cohort performance do not generalize in the cross-cohort context i.e., 347 cohort specific predictors are being learned in the within-cohort context. 348 349 The five overlapping genes between the within and cross-cohort approaches were RASSF6, BMP1, MXRA8, 350 SAMD11 and LLGL1. RASSF6 encodes a Ras-effector protein which binds activated Ras and is involved in 351 induction of apoptosis . It’s a member of a tumor repressor network which is highly conserved , and it’s 352 inactivation or suppression has been linked to several diseases including cancer (22). BMP1 (bone 353 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint morphogenic protein 1) is a metalloprotease which regulates formation of extracellular matrix (ECM) by 354 cleaving/processing of various precursor proteins. Mutations in this gene have been associated with 355 skeletal/connective tissue disorders (23). The role of MXRA8 (Matrix Remodeling Associated 8- encoding 356 gene) is less established, but it is thought to be involved in glial blood -brain barrier functioning and has 357 been associated with glioma progression and a suggested therapeutic target for glioma (24). MXRA8 is 358 known to play roles in inflammation and cell signaling that could be potentially relevant to endometriosis. 359 Long non-coding RNA transcribed from the SAMD11 region has been functionally linked to endometrial 360 biology (decidualization) (25). LLGL1 is a well -known tumor suppressor gene involved in maintaining 361 epithelial cell polarity and tissue architecture in various cancers, including endometrial cancer (26). 362 363 Classifier-aligned integrated gradients (CA -IG) allows FM embedding gene–level 364 interpretation 365 366 Next, we sought to understand if and how model explainability differed when using FM embeddings for 367 disease classification. Because the downstream classifier is not jointly fine-tuned with the encoder, standard 368 gradient attribution cannot be directly applied. We therefore developed classifier -aligned integrated 369 gradients (CA -IG), an attribution approach that assigns importance to FM gene tokens by computing 370 Integrated Gradients (IG) through the frozen encoder along a classifier-relevant direction in the embedding 371 space. This direction is defined locally for each sample using SHAP scores from the downstream classifier 372 with respect to the FM embedding, ensuring that gene -level attributions reflect the features used by the 373 classifier. Importantly, SHAP is used solely to define this classified -aligned weighting of embedding 374 dimensions, while IG performs the attribution itself. This formulation avoids backpropagation through the 375 non-differentiable classifier, substantially reducing computational cost while preserving biological 376 interpretability by providing gene-level explanations (see Methods). 377 378 379 Table 5. FM approach. Top 20 explainable genes identified via CA-IG for Geneformer embedding model 380 for the within-cohort and cross-cohort approach. SHAP scores normalized using a positive-only approach, 381 then averaged across runs and folds (HGNC symbols). Asterisks (*) represent genes found by both 382 approaches. (Full gene list available in Supplementary Tables S9 (within-cohort) and S10 (cross-cohort)). 383 384 WITHIN COHORT CROSS-COHORT GeneID CA-IG score Gene name GeneID CA-IG score Gene name .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint DDIT3* 0.01254973 DNA Damage Inducible Transcript 3 DDIT3* 0.00966168 DNA Damage Inducible Transcript 3 LRRC3C* 0.00604412 leucine rich repeat containing 3C LRRC3C* 0.0058623 leucine rich repeat containing 3C TBC1D3F* 0.00598382 TBC1 domain family member 3F TBC1D3F* 0.00584974 TBC1 domain family member 3F OR1J2* 0.00528166 olfactory receptor family 1 subfamily J member 2 OR1J2* 0.00370744 olfactory receptor family 1 subfamily J member 2 FRG2* 0.00337272 FSHD region gene 2 FRG2* 0.00330829 FSHD region gene 2 TBC1D3E* 0.00240315 TBC1 domain family member 3E MEPCE* 0.00101443 methylphosphate capping enzyme PPIAL4F* 0.0015164 peptidylprolyl isomerase A like 4F TBC1D3E* 0.00099277 TBC1 domain family member 3E MEPCE* 0.00144943 methylphosphate capping enzyme TBC1D3K* 0.0009618 TBC1 domain family member 3K OR51B5 0.00127887 olfactory receptor family 51 subfamily B member 5 TBC1D3* 0.00094298 TBC1 domain family member 3 TBC1D3K* 0.00124146 TBC1 domain family member 3K RSC1A1* 0.00092813 regulator of solute carriers 1 TBC1D3C* 0.00122356 TBC1 domain family member 3C PPIAL4F* 0.00089612 peptidylprolyl isomerase A like 4F RABGAP1* 0.00120914 RAB GTPase activating protein 1 RABGAP1* 0.00089356 RAB GTPase activating protein 1 RSC1A1* 0.00119487 regulator of solute carriers 1 DERPC 0.00088404 DERPC proline and glycine rich nuclear protein TBC1D3* 0.00118897 TBC1 domain family member 3 SPDYE14* 0.00075791 speedy/RINGO cell cycle regulator family member E14 TBC1D3H* 0.00116252 TBC1 domain family member 3H TBC1D3C* 0.00074984 TBC1 domain family member 3C SPDYE10* 0.00109991 speedy/RINGO cell cycle regulator family member E10 PCDHGC5* 0.00073012 protocadherin gamma subfamily C, 5 SPDYE15 0.00100317 speedy/RINGO cell cycle regulator family member E15 TBC1D3H* 0.0007123 TBC1 domain family member 3H PCDHGC5* 0.00097082 protocadherin gamma subfamily C, 5 PIGY 0.00069504 phosphatidylinositol glycan anchor biosynthesis class Y PARP2* 0.00096208 poly(ADP-ribose) polymerase 2 PARP2* 0.0006938 poly(ADP-ribose) polymerase 2 SPDYE14* 0.00093309 speedy/RINGO cell cycle regulator family member E14 SPDYE10* 0.00068962 speedy/RINGO cell cycle regulator family member E10 385 386 387 Using our CA -IG approach to extract explainability for the FM -embedding classifiers, we observed a 388 striking conservation of disease -associated genes between the within -cohort and cross -cohort settings 389 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint (Table 5). Of the top 20 genes identified in the within -cohort analysis, 18 overlapped within the cross -390 cohort top 20, and conversely 18/20 of the top cross-cohort genes were also present in the within-cohort list 391 (Table 5). Notably, the top five genes were identical and ranked in the same order for both strategies, 392 indicating highly consistent relative importance across evaluation regimes (Fig 4). This contrasts sharply 393 with the TPM-based models, where predictive features varied substantially between within and cross cohort 394 analyses. The stability of the embedding -derived attributions parallels the reduced performance drop 395 observed for the FM-based classifiers under cohort shift, suggesting that the pretrained representations may 396 prioritise transferable disease-related signal over cohort-specific effects. 397 398 The five most predictive genes across both settings were DDIT3 (DNA Damage-inducible transcript 3 -399 encoding gene), LRRC3C (Leucine Rich Repeat Containing 3C-encoding gene), TBC1D3F (TBC1 Domain 400 Family Member 3F), OR1J2 (Olfactory Receptor Family 1 Subfamily J Member 2 -encoding protein) and 401 FRG2 (FSHD Region Gene 2-encoding protein) (Fig 4). 402 403 404 405 406 Figure 4: Conserved gene-level attributions from Geneformer embeddings. Ten highest–scoring genes 407 for the frozen Geneformer encoder, computed with Integrated Gradients (PAD baseline; Methods) and 408 aggregated over true-positive test samples, runs/seeds, and outer folds. Panels show the within-cohort (left) 409 and cross-cohort (right) strategies. Values are per -sample positive-only, L1-normalised attribution scores 410 (positive SHAP values kept only), with larger indicates more disease-driving evidence. 411 WITHIN-COHORT CROSS-COHORT .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 412 The most important gene for endometriosis disease prediction using both the within-cohort and cross-cohort 413 approaches was DDIT3, which encodes a transcription factor which mediates response to cellular stress , 414 with a main role in promoting apoptosis during conditions of severe stress. Endoplasmic reticulum (ER) 415 stress is being increasingly recognized as playing a major role in the pathogenesis of endometriosis (27,28), 416 and DDIT3/CHOP is a crucial transcription factor activated by ER stress which has also been shown to 417 further upregulate ER stress and apoptosis in endometriosis disease (29). DDIT3 promotes apoptosis (30), 418 which may drive the degeneration of ectopic tissue. DDIT3 has been strongly linked to a variety of other 419 diseases including inflammatory diseases, cancer, neurodegenerative diseases and metabolic diseases 420 (31,32). Its’ overexpression has been correlated with poor breast cancer prognosis and implicated in disease 421 progression via its role in cellular stress response pathways (33). DDIT3 plays key roles in inflammatory 422 diseases by activation of several inflammatory pathways including NF-kb pathway (34) and is upregulated 423 by cytokines TNF -a and IL-1B. DDIT3 overexpression has also been linked to reduced antiviral innate 424 immune activity, encouraging viral replication and disease (35). This gene is a strong candidate for further 425 exploration and validation. 426 427 The second ranked gene, LRRC3C, encodes a (predicted) membrane -located protein containing leucine -428 rich repeats motifs (LRR), important for protein-protein interactions. Although mechanistic studies linking 429 LRRC3C with disease remain limited, LRRC3C has been implicated as important in inflammatory disease, 430 particularly in inflammatory bowel disease susceptibility via its impact on cell proliferation and apoptosis 431 (36). 432 433 The protein encoded by the third most important gene TBC1D3F is one of eight paralogs (A-H) belonging 434 the TBC1D3 family. Members of this family of molecules have been found to be over expressed in several 435 cancer types (prostate, breast, pancreatic) (37) and have been associated wi th poor prognosis of kidney 436 renal cell carcinoma (38). TBC1D3F was found to be overexpressed in the epithelial ovarian cancer cell 437 lines and primary tumors compared to normal tissues (39). Alongside TBC1D3F, several members of the 438 TBC1D3 family are present in the top 20 list of most important genes for endometriosis disease prediction, 439 including TBC1D3E, TBC1D3K and TBC1D3RC. The TBC1D3 locus is frequently associated with cancer, 440 amplified in 15% of prostate cancers , and the TBC1D3 cluster is extremely structurally variable (40,41). 441 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint The amplified presence of genes from the TBC1D3 family ( 6/20 of the topmost important genes for both 442 within and cross-cohort models) makes this family a strong candidate for experimental validation. 443 444 Olfactory receptors (ORs) are expressed in several non -olfactory tissues including the endometrium and 445 are linked to inflammation through specific signaling pathways. Although information on OR1J2 is 446 relatively sparse, other ORs, for example, OR5H2, have been identified as a potential target in endometrial 447 cancer biology related to the IGF1 signaling pathway, and may be involved in cancer protection in some 448 individuals (42). 449 450 LLM-driven knowledge graph inference reveals enriched biological processes linked to 451 endometriosis 452 453 Having first investigated the individual functional roles of the top genes and their potential relevance 454 endometriosis disease pathology, we next sought to characterise the shared molecular landscape of these 455 genes to determine whether they shared any links or associations with each other and/or with endometriosis. 456 Because the molecular mechanisms of endometriosis remain poorly characterised, we employed a 457 literature-derived knowledge-graph and LLM-inference approach to do this, which may reveal relationships 458 which traditional bioinformatic gene-centric analysis may miss. For each of the five most important genes 459 for endometriosis prediction (DDIT3, LRRC3C, TBC1D3F, OR1J2 and FRG2), we systematically retrieved 460 all PubMed-indexed articles in which it was mentioned and identified semantic relationships between that 461 gene and identified biomedical entities such as disease, gene, pathway, proteins etc. All entities and 462 relationships were represented as entity -entity-relation triples and integrated into a unified bipartite 463 knowledge graph. 464 465 The resulting knowledge graph contained 312,598 triples linking 81,219 distinct entities of which 9,835 466 were proteins, 8,669 were diseases, 5,800 were genes, 4,698 were biological process and 1,367 were drugs 467 (Full list entity types in Supplementary Table S11). The graph contains 33,092 distinct relation types with 468 the most common being co-expression (4,910 triples) and regulation (4,095 triples) (Full list relation types 469 in Supplementary Table S12). First, we investigated the similarity of connected entities and relations 470 between the five important genes by calculating i) normalized shared -entity overlap, and ii) relation-style 471 similarity for each gene pair (see Methods). As expected, most gene pairs exhibited relatively low 472 normalized shared-entity overlap, suggesting genes tended to engage with largely non -overlapping entity 473 sets (Supplementary Table S13). Due to the vastly larger number of triples extracted for gene DDIT3 than 474 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint the other important genes, 27%-43% of unique entities from genes LRRC3C, TBC1D3F, OR1J2 and FRG2 475 were shared with DDIT3 (Supplementary Table S13). 476 477 To investigate how our five important genes may converge on common biological mechanisms linked to 478 endometriosis disease pathology, we reduced the multi-gene knowledge graph to a compact JSON structure 479 containing shared entities, shared relations and inter -gene similarity score, and provided as structured 480 context to an LLM to suggest biologically plausible mechanisms supported by the evidence. Although the 481 five genes have distinct primary functions, the LLM results highlighted their roles in several convergent 482 biological themes linked to endometriosis. ER stress is being increasingly recognized as playing a major 483 role in the pathogenesis of endometriosis (27,28), and is enhanced by DDIT3/CHOP in endometriosis (29). 484 The LLM results correctly characterize DDIT3s potential role in endometriosis disease and links it to the 485 role of TBC1D3F in vesicle-mediated cytokine release (43), suggesting a shared influence on inflammatory 486 and stress -response pathways in endometriosis. LLM results also suggest a role of TBC1D3F in 487 endometriosis pathogenesis via its role in dysregulated vesicle trafficking, which can alter secretion of 488 proinflammatory cytokines and growth factors that promote ectopic endometrial implantation and 489 contribute towards the pro -inflammatory milieu. This is supported by the literature which indicate that 490 extracellular vesicles play a role in endometriosis disease pathogenesis (44). LLM results suggest a joint 491 contribution of genes LLRC3C and FRG2 to endometriosis pathogenesis through processes related to 492 enhanced adhesion and fibrosis, mediated by extracellular -matrix remodeling ( LLRC3C) and fibroblast 493 proliferation/fibrosis (FRG2). In the literature LRRC3C is only partially characterized with predicted roles, 494 however LRR (leucine-rich repeat) proteins have been implicated in cell-adhesion and extra-cellular-matrix 495 (ECM) remodeling (45,46). Enhanced adhesion of endometrial stomal cells to peritoneal surfaces and 496 altered ECM are key steps in endometrial lesion establishment (47,48). 497 498 BMFM attribution rankings for GSEA reveals enriched gene pathways linked to 499 endometriosis pathophysiology 500 501 To further investigate associations between the model -derived attribution scores and potential 502 endometriosis disease mechanisms, we performed gene set enrichment analysis (GSEA) from the cross -503 cohort embedding descending gene list ranked by CA -IG score (see Methods). This analysis identified 504 several positively enriched pathways (FDR < 0.05), indicating a significant association with disease 505 classification. Enriched disease pathways were predominantly related to inflammatory cytokine signaling, 506 cell stress and cell survival, mechanisms consistent with the current understanding of endometriosis disease 507 pathophysiology (Supplementary Figure 1, Supplementary Tables S14 -S17). Within cytokine signaling, 508 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint pathways responsible for positive regulation of type II interferon (IFN-gamma), interleukin-17 (IL-17) and 509 interleukin-8 (IL-8) were all significantly enriched amongst top ranking genes (FDR < 0.05, Supplementary 510 Figure 1). Recent evidence has suggested IL -17 may be an important mediator of endometriosis disease, 511 with elevated levels observed across several endometriosis studies (49,50) and sample types including 512 endometriotic lesions (51) and peritoneal fluid (52) from women with endometriosis. Furthermore, higher 513 levels of IL -17-producing cells have been found in the peripheral blood of women with endometriosis, 514 accompanied by increased IL -17 protein expression in their endometriotic tissue (53), indicating both 515 systemic and local involvement of IL -17 in endometriosis disease. IL -8 is also frequently elevated in the 516 peritoneal fluid of women with endometriosis, and in vitro studies have demonstrated IL-17 stimulation of 517 IL-8 secretion from endometriotic stomal cells, linking both pro -inflammatory cytokines to sustenance of 518 the pro-inflammatory milieu of endometriosis (54). Elevated levels of type II interferon (IFN-gamma) have 519 also been detected in endometriotic tissue and peritoneal fluid of women with endometriosis, supporting a 520 role for IFN -gamma mediated inflammatory and immune dysregulation in endometriosis disease 521 pathophysiology (55,56), however some studies have reported conflicting results (57). Alongside cytokine 522 signaling, GSEA indicated several enriched pathways associated with disease were linked to cell stress and 523 cell survival mechanisms, including endoplasmic reticulum (ER) stress, unfolded protein response (UPR) 524 and regulation of apoptosis signaling pathways, highlighting altered cell stress and survival programs in 525 endometriosis (Supplementary Figure 1). ER stress and UPR activation have been shown in endometriotic 526 tissue, and is thought to reflect a response of ectopic cells to the pro -inflammatory milieu (28,58,59). 527 Enrichment of pathways regulating ER stress-induced apoptosis is consistent with a dysregulated balance 528 between stress responses and programmed cell death, a process which has been widely implicated in 529 endometriotic lesion persistence via survival of endometrial cells via evasion of apoptosis (60,61). Together 530 with existing literature, these findings support the interpretation that stress-altered regulation of cell survival 531 pathways represent a prominent feature of endometriosis pathophysiology. 532 533

534 535 In this study, we leveraged a multi-cohort bulk RNA -seq dataset to systematically assess the benefit of 536 RNA foundation model (FM) embeddings as feature representations for endometriosis disease prediction 537 using a downstream classifier. We evaluated the benefits of these FM-derived features using a cohort aware 538 strategy and found that FM-derived representations consistently and substantially outperformed gene-level 539 TPM features under a cross -cohort evaluation approach, where models were required to generalize to 540 entirely unseen datasets. This setting reflects a common challenge in clinical transcriptomics and biomarker 541 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint discovery, and our results demonstrate that pre-trained RNA representations improve robustness to cohort 542 heterogeneity. 543 544 A central contribution of this work is the introduction of classifier-aligned integrated gradients (CA-IG), a 545 computationally efficient approach for interpreting gene-derived FM embeddings. While IG benefits from 546 a GPU for speed, it substantially lowers the computational burden for interpretability of FMs and makes 547 analyses feasible with modest resources. Applying CA -IG to Geneformer embeddings revealed highly 548 stable gene-level attributions across both cohort-aware strategies, in contrast to TPM-based models where 549 predictive genes varied substantially between cross and within cohort approaches. Embedding models 550 consistently prioritised a conserved set of genes, supporting the notion that pretrained encoders can capture 551 generalisable biological knowledge rather than cohort -specific noise. Pathway -level analyses further 552 revealed convergence on inflammatory cytokine signaling, cell stress and cell survival mechanisms, 553 processes strongly implicated in endometriosis disease pathophysiology. 554 555 Some limitations remain within this work. First is the small size of several of the individual cohorts, which 556 increased variance in train/test splits. We addressed this using dataset-aware splits, a grouped-cohort outer 557 CV, a nested inner CV for model selection, and multiple seeds with bootstrap confidence intervals. Second 558 was the imbalanced nature of the full cohort, with comparatively few controls. This imbalance likely 559 contributed to the reduced predictive performance of class 0 compared to class 1 and is a common issue in 560 disease/control classification . Our third limitation was associated with inconsistent reporting of key 561 covariates due to small and varied sample sizes (e.g., menstrual -cycle phase, lesion subtype), which 562 prevented adjustment for potential confounders. 563 564 RNA FMs offer an accessible and powerful solution for improving the robustness of disease classification 565 in transcriptomics. Here we demonstrated this for endometriosis, a condition hampered by diagnostic delays 566 and biological heterogeneity, where this CPU-friendly approach provided more generalizable predictions 567 and stable, interpretable biomarkers, laying a practical foundation for future application to other diseases 568 and use-cases, leveraging large-scale biological priors in clinically focused research. 569 570 571

572 573 Datasets and pre-processing 574 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 575 We assembled a multi -cohort bulk RNA-seq corpus from the Gene Expression Omnibus (GEO) (20) by 576 searching for human studies of endometrial tissue related to endometriosis and restricting results to 577 expression profiling by high-throughput sequencing. Single-cell, spatial transcriptomic studies, and studies 578 with fewer than 10 samples (pre -filtering) were excluded. We used the GEO identifiers for our selected 579 studies to download the associated transcriptomic data for these studies and firstly, to merge and normalise 580 the gene expression count data (TPM) into a unified data matrix, and then to extract harmonised metadata 581 for the same samples i.e., in this case we used an LLM to yield a label for each sample of endometriosis 582 disease yes (class 1) or control/no (class 0). Labels were derived from the harmonised metadata and then 583 further checked and curated by an expert molecular biologist to remove ambiguous cases (e.g., endometrial 584 cancer, unclear disease status). After filtering, the final corpus comprised 12 studies and 334 samples (259 585 endometriosis, 75 healthy controls). GEO Series IDs were retained as a grouping variable for subsequent 586 cohort-aware cross-validation. Raw counts were converted to transcripts per million (TPM) using GRCh38 587 (p13) gene lengths and log-transformed (log(1+TPM)). Gene identifiers were harmonised across cohorts, 588 and all matrices were aligned to a common 17,518-gene universe. 589 590 Fold construction 591 592 Within-cohort 593 594 For the within -cohort analysis, we used a dataset -aware, stratified 3 -fold split approach. To ensure each 595 dataset was represented in the training set for each fold, we first selected one sample per GEO series or 596 dataset as the ‘anchor’ sample and pinned all anchors to the training set in every fold. We then applied 597 stratified-K-fold (K=3, shuffled, fixed seed) to the non -anchor samples only, preserving the global class 598 ratio in each test fold. In fold k the training set was the union of all anchors and the non-anchor samples not 599 assigned to the k-th test split; the test set was the corresponding stratified k-th split. This procedure yields 600 zero sample overlap between train and test, keeps test-set class balance approximately intact, and enforces 601 an in-domain condition in which any series or dataset appearing in test also has at least one sample (its 602 anchor) in train. Fold assignments and per-fold counts are reported in Table 6. 603 604 Cross-cohort 605 606 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint For the cross-cohort approach, we constructed balanced grouped-cohort splits with K=3 outer folds to match 607 the within-cohort approach. Entire datasets (GEO Series IDs) were assigned to test folds by a constrained 608 random search (1,000 trials; fixed seed) that (i) enforced a minimum test size (at least 60 samples per fold) 609 and (ii) minimised the across-fold standard deviation of the positive fraction and of fold sizes. This resulted 610 in test/train dataset assignments as reported in Table 6. 611 612 613 Table 6. Fold construction for the within-cohort and cross-cohort folds. 614 615 Within- cohort Folds Train (n) Test (n) Positive proportion (train) Positive proportion (test) Dataset IDs 1 226 108 0.78 0.77 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 2 227 107 0.78 0.77 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12 3 227 107 0.78 0.78 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12 Cross- cohort Folds Train (n) Test (n) Positive proportion (train) Positive proportion (test) Dataset IDs 1 237 97 0.78 0.77 6,7,9 2 228 106 0.77 0.78 5,8,12 3 203 131 0.78 0.77 2,4,10,11 616 617 Extraction of Foundation Model Embeddings 618 619 We evaluated five pretrained encoders spanning bulk and single -cell transcriptomic pretraining: 620 BulkRNABERT (18), scFoundation (17), scGPT (16), Geneformer (15) and BMFM-MLM-RDA (19). 621 These models learn transferable structure from large expression corpora and produce compact, context -622 aware representations (embeddings) for downstream prediction. 623 624 For each sample we performed feature harmonisation as follows: (i) mapping the TPM vector to the 625 encoder’s gene vocabulary (e.g. HGNC to Ensembl); (ii) dropping genes absent from that vocabulary; (iii) 626 consolidating duplicate symbols by averaging their TPM values; and (iv) constructing model -specific 627 feature sets (tokens). To extract embeddings per sample, we then executed a single forward pass with the 628 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint encoder in evaluation mode (frozen weights) and obtained a fixed -dimensional embedding per sample. 629 Embeddings were computed once per sample and cached before any model selection; identical splits and 630 seeds were reused across models for paired comparisons. We intentionally avoided additional downstream 631 scaling or feature selection on embeddings to preserve the pretrained geometry. When multiple HGNC 632 symbols mapped to the same Ensembl identifier, we averaged their TPM values by grouping them by their 633 shared Ensembl identifier prior to input construction. We chose this approach as averaging preserves scale 634 and avoids artificial inflation for genes with many aliases, which could otherwise bias downstream models. 635 Model-specific approaches are detailed below: 636 637 Geneformer 638 Geneformer is a transformer pretrained on 104M single-cell transcriptomes with rank–value encodings and 639 masked-token objectives (15). We mapped TPM(log+1) matrices to Ensembl IDs, tokenised with the 640 official toolkit, and extracted the penultimate layer. Mean pooling excluded special tokens (CLS/EOS, 641 PAD) producing a 768-dimension embedding per sample. 642 643 BulkRNABERT 644 This encoder-only model was pretrained on bulk RNA -seq (GTEx/ENCODE/TCGA) with a masked-LM 645

over discrete expression tokens and permutation-invariant gene order (18). TPM(log+1) matrices 646 were harmonised to the accepted Ensembl set and ordered to the model list. We extracted embeddings from 647 the final layer mean-pooled across tokens to yield a 256-dimensional embedding per sample. 648 649 scFoundation 650 ScFoundation is an asymmetric encoder–decoder pretrained on >50 million human single-cell profiles (62). 651 TPM(log+1) vectors were aligned to the 19,264 -gene vocabulary (missing genes padded with zeros). The 652 official routine returns a 3,072 -dimensional vector comprising four 768 -dimensional pools (1: last token, 653 2: penultimate token, 3: max pooled, 4: mean pooled). We split this file and used the last 768 dimensions 654 for all further analysis (mean-pooled 768-dim block). 655 656 scGPT 657 ScGPT is a generative transformer pretrained on 33 million single -cell profiles with masked -attention 658 objectives, jointly learning gene and cell embeddings (16). TPM(log+1) matrices were aligned to the model 659 vocabulary and embedded, yielding a 512-dim vector per sample. 660 661 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint BMFM-MLM-RDA 662 We used the MLM -RDA variant of BMFM: a BERT -base encoder, pretrained with read -depth–aware 663 (RDA) down sampling pretrained via expression masking and RDA downsampling 664 (https://huggingface.co/ibm-research/biomed.rna.bert.110m.mlm.rda.v1, (19)). We used TPM (log+1) data 665 for a single forward pass through the frozen embedder; the final -layer [CLS] vector (768 -dimensional) 666 served as the sample embedding. 667 668 Explainable Machine Learning workflow 669 670 Train and test protocol 671 We used the open-source explainable-AI tool AutoXAI4Omics (https://github.com/IBM/AutoXAI4Omics) 672 to run a robust classical explainable ML analysis on both TPM and embedding inputs for both the within -673 cohort and cross-cohort strategies, ensuring consistency, transparency, and reproducibility across model 674 comparisons (63). We chose the AdaBoostClassifier model as our downstream light-weight classifier (64). 675 676 We evaluated two input data types within the same AutoXAI4Omics classical -ML pipeline for 677 endometriosis classification: (i) log TPM values and (ii) foundation-model (FM) embeddings derived from 678 log TPM values. We adopted two strategies; i) within-cohort, where train and test contain samples from the 679 same GEO dataset ID, and ii) cross -cohort, where the held -out test fold comprised unseen datasets not 680 represented within the training set. The outer loop used K=3 folds; in each outer fold, models were trained 681 on K-1 folds and evaluated once on the held-out fold. Within the fold's training set, hyperparameters were 682 tuned by stratified 5-fold cross-validation, selecting the configuration that maximised weighted-F1 (primary 683 metric), then refitting that configuration on the full training data before testing. To assess robustness, each 684 outer split was repeated with five different random seeds; we report means across folds and seeds. All 685 preprocessing (scaling, feature transforms, encoders and feature selection) was fitted strictly on the training 686 side and applied unchanged to validation/test to avoid leakage. 687 688 TPM data preparation for ML 689 For all AutoXAI4Omics runs with TPM data we performed feature selection using the following two inbuilt 690 methods; (i) first a variance filter removed near -constant genes (threshold = 1); and (ii) secondly a 691 univariate ANOVA test (F-statistics) was used to select the top k genes, where k was chosen automatically 692 from a base-10 grid with a minimum of 25 features (no explicit maximum). For each candidate k, a small 693 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint random forest was trained on the selected genes and evaluated by F1 score. From the set of values yielding 694 high performance, we selected the smallest k, prioritising a smaller feature set while retaining performance. 695 696 FM-embedding data preparation. 697 For all FM vectors we disabled standardisation and feature selection, treating embeddings as fixed 698 representations to avoid distorting pretrained geometry. 699 700 Post-AutoXAI4Omics Bootstrapping and significance calculation 701 702 To obtain a robust and sound estimate of all models performance variability and allow statistically valid 703 comparison of models, we calculated confidence intervals for fold -averaged metrics via hierarchical 704 bootstrapping, with 10,000 iterations. For each model we had performance metrics from 3 folds and 5 seed 705 runs per fold, totaling 15 runs per model. Each bootstrap iteration consisted of drawing 3 folds with 706 replacement, and per selected folds drawing 5 seed runs with replacement. For each metric ‘fold means’ 707 were calculated which were then used to calculate ‘overall means’ across all folds. The same approach was 708 used to calculate the effect size (while preserving paired seeds across models), defined as the difference in 709 model metrics to the baseline model (TPM input, not embeddings). Bootstrap estimates of the mean and 710 effect size over 10,000 iterations were used to construct 95% confidence intervals. To test whether 711 performance of paired models differed significantly we leveraged the same nested -bootstrap samples that 712 we used for confidence intervals. In each of the 10,000 bootstrap iterations we paired the corresponding 713 runs of two models (same index), computed the difference in performance metric for each fold, and 714 averaged these differences over the 3 folds. This yielded one bootstrap difference value per metric per 715 iteration. Under the null hypothesis that the models have equal performance, the distribution of differences 716 would be centred on zero (no difference). We therefore we calculate two -tailed p-values as follows: i) 717 compute the proportion of bootstrap differences which lie above zero, and the proportion that lie below 718 zero, ii) take the smaller proportion (representing the least likely direction of difference) and iii) we double 719 this to get a 2 -tailed p-value. This approach was applied to all pairwise comparisons (baseline/vanilla vs. 720 all BMFMs) for the metrics F1-score, F1-score (class 0), F1-score (class 1) and accuracy. 721 722 Evaluation Metrics 723 724 Primary reported metrics were weighted -F1, F1 class0, F1 class1 and accuracy with 95% CIs obtained via 725 bootstrapping as described above. To characterize the benefit of frozen FM embeddings over the baseline 726 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint TPM we report paired effect sizes for weighted-F1, F1class0, F1class1 and accuracy, computed paired by seed 727 within fold and summarized with the same bootstrap (paired 95% CIs). 728 729 Explainability 730 731 SHAP-based gene attribution (TPM models). 732 733 We used Shapley Additive exPlanations (SHAP) [ 26, 27] as implemented in AutoXAI4Omics to quantify 734 per-gene contributions to the positive-class F1 score. We used the TreeSHAP explainer for the AdaBoost 735 classifier, fitted on training folds only and applied to the held -out test split. To focus on disease -driving 736 signal (rather than model error or protective evidence), we collected and aggregated SHAP scores from all 737 folds (n=3) and all runs (n=5) of the AdaBoost classifier only for true-positive test samples (ground-truth 738 is ‘disease’ and predicted positive threshold >= 0.5). To allow comparison of SHAP scores across samples, 739 aggregated true-positive SHAP scores were subject to a positive-only normalization, where positive SHAP 740 scores were normalized per sample by scaling each sample’s positive SHAP vector into a distribution 741 summed to 1(only positive SHAP scores were retained to allow focus on features driving disease 742 prediction) Normalized SHAP values were then averaged across folds and seeds to generate a normalized 743 average SHAP value per feature. 744 745 Geneformer attribution (classifier-aligned integrated gradients) 746 747 Focusing on our best performing FM Geneformer, in this study we used transcriptomic information as input 748 into the FM to derive embeddings (see Methods: Extraction of FM embeddings, Geneformer) . We then 749 used these embeddings to train a downstream classifier (AdaBoost) to predict endometriosis status , but 750 notably this classifier was not fine-tuned jointly with Geneformer. Using explainability techniques, our goal 751 was to trace disease predictions back through the FM encoder to identify which input genes were most 752 responsible for the classification decision from our resultant AdaBoost model. To derive this explainability, 753 we developed an approach that we call classifier-aligned integrated gradients (CA-IG) to attribute the FM 754 token inputs (gene tokens) to the downstream disease score used for classification using Integrated 755 Gradients (IG), a well-established attribution method (65). However, the “classifier-aligned” aspect of our 756 approach, computing gradients using SHAP values rather than through the entire model, brings novelty and 757 utility, dramatically reducing computational cost while maintaining gene -level information about what 758 drives predictions. 759 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 760 When a tokeni sed RNA-seq sample is processed by Geneformer , it produces a matrix of dimension 761 𝐿 × 𝐷_𝑚 (where L represents length and D_m is the embedding size). The Geneformer encoder at layer ℓ 762 (the penultimate layer in all our experiments) transforms this into hidden states of dimension L×D. We 763 focused on real tokens only, excluding special tokens like CLS, EOS, and PAD, and appl ied a masking 764 procedure to obtain the sample embedding. Therefore, the sample embedding, denoted as 𝑧(𝑋), was created 765 by mean-pooling the masked hidden states over all real tokens, weighted by their mask values 𝑚 ∈ {0,1}^𝐿. 766 This produced a D-dimensional vector (768 dimensions in this instance) that represented the entire sample 767 in the embedding space. 768 769 Our downstream classifier AdaBoost operates on this D -dimensional embedding space to produce a 770 positive-class score (the disease probability). Since the classifier is not fine-tuned jointly with Geneformer 771 we used a sample-specific surrogate direction derived from the reali sed SHAP attribution vector of the 772 classifier's positive logit with respect to this embedding. This surrogate direction 𝑤(𝑥) exists in the same 773 D-dimensional space as the embedding. We then defined a scalar objective 𝑠(𝑋) as the inner product 774 between this SHAP direction and the sample embedding 𝑧(𝑋). This scalar objective is locally faithful to 775 the classifier and enable d efficient backpropagation into the token embeddings without requiring 776 differentiation through the separate classifier. We call this approach classifier-aligned integrated gradients 777 (CA-IG) because it aligns the gradient direction with the classifier's behaviour. 778 779 Computationally, this approach converts the vector sensitivity into a single directional derivative along the 780 SHAP direction 𝑤(𝑥). Using the linearity of gradients, computing the gradient of this scalar objective with 781 respect to the embedding X can be expressed as a weighted sum of per -dimension gradients, where each 782 dimension is weighted by the corresponding SHAP weight. This means a single backward pass per IG step 783 yields the same weighted gradient that a user would obtain by computing per -dimension gradients 784 separately and weighting them afterward. This reduces the computational cost from 𝑂(𝑆 · 𝐷) to 𝑂(𝑆) 785 backpropagations for S integration steps across embedding dimension D. To prevent data leakage, the 786 SHAP direction 𝑤(𝑥) is computed using only the training fold, then applied to held-out test samples. 787 788 IG is a path-based attribution method that accumulates gradients along a straight-line path from a baseline 789 to the actual input. We used a PAD baseline B of dimension 𝐿 × 𝐷_𝑚, where each position was set to the 790 learned pad token vector. The IG from baseline B to input X were calculated by integrating the partial 791 derivative of the objective with respect to the input along the path from B to X, parameterised by α from 0 792 to 1. This continuous integral was approximated using a Riemann sum with S=32 discrete steps along the 793 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint path. At each step, we evaluated the gradient at a point that interpolated between the baseline and the actual 794 input. Token-level scores were obtained by summing these IG across all channels (the D_m embedding 795 dimensions) for each token position, resulting in a score 𝑎_𝑖^ + (𝑥) for each gene token i. These token 796 scores were then mapped to genes by removing special tokens and using the tokeniser's gene dictionary to 797 form gene-level attributions 𝑔^ ± (𝑥). 798 799 Akin to the SHAP analysis, we focused on disease-driving signals by reporting summaries that emphasised 800 positive (disease-promoting) contributions. For each sample x, we normali sed the attributions per sample 801 before aggregating. We created a positive-only summary by taking the maximum of each gene's attribution 802 and zero, then normali sing these positive values to sum to 1 (with a small constant ε added to the 803 denominator for numerical stability). This g ave us ĝ_𝑖^ + (𝑥), representing gene i's relative positive 804 contribution to the disease prediction for sample x. We also compute d a signed summary ĝ_𝑖(𝑥) that 805 preserved the direction of effects by normali sing using the sum of absolute values across all genes, again 806 with the stabilising constant ε. We attributed only true-positive patients—specifically, test samples where 807 the ground truth y=1 (disease present) and the predicted probability exceeded 0.5. For outer CV fold k and 808 repeat r, we identif ied the set 𝑇_𝑘, 𝑟 of true -positive test samples. We then average d the normali sed 809 attributions within each run, within each fold, and finally across all folds. 810 811 The final importance score 𝐼𝐺𝐼𝑚𝑝_𝑖^ + for gene i was calculated by averaging over K folds, R repeats, and 812 all true-positive samples within each fold -repeat combination, using the positive -only normalised scores. 813 Similarly, 𝐼𝐺𝐼𝑚𝑝_𝑖^ ± used the signed normalised scores to capture directional effects. This multi-level 814 averaging produced stable, generali sable gene importance scores that reflect consistent patterns across 815 different data splits and random initialisations, minimising the influence of any single model or sample. To 816 our knowledge, CA -IG is the first attribution protocol tailored to frozen RNA foundation -model 817 embeddings with a separate downstream classifier, and it reduces the per-sample IG complexity from O(S 818 ·D) to O(S) while remaining CPU-efficient. We report IGImp+ in the main text and provide IGImp± in the 819 Supplementary material. 820 821 GSEA from attribution ranks 822 823 To identify enriched biological pathways associated with disease we performed gene set enrichment 824 analysis (GSEA) using pre-ranked gene lists derived from feature importance scores. Four separate analyses 825 were performed with results from the following approaches: a) within-cohort TPM model and associated 826 gene SHAP scores; b) cross -cohort TPM model and associated gene SHAP scores; c) within-cohort FM 827 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint embedding model (Geneformer) and associated gene attribution scores from CA -IG and d) cross -cohort 828 FM embedding model (Geneformer) and associated gene attribution scores from CA -IG. For TPM-based 829 analyses gene importance scores were obtained from SHAP values of the disease classification model whilst 830 FM GSEA analyses were run with gene attribution scores calculated from the classifier-aligned integrated 831 gradients method. The gene universe was defined as all genes present in the input TPM matrix. For each of 832 the four analyses, gene importance scores were mapped to the gene universe and genes were ranks in 833 descending order of their signed importance. Pre-ranked GSEA was performed using GSEApy (66) with 834 10,000 permutations, testing using gene se ts from the Enrichr GO Biological Processes library (67). 835 Pathways significantly enriched at either end of the ranked gene list were identified at FDR < 0.0 5, with 836 positive normalized enrichment score (NES) indicating enrichment among disease -associated genes and 837 negative normalized enrichment score (NES) indicating enrichment among healthy -associated genes. 838 Parameters used are as follows : min_size= 15, max_size= 500, permutation_num= 10,000, Benjamini –839 Hochberg FDR control (seed= 42). For each term we report the NES, nominal P value, and FDR-corrected 840 P-value. “Leading-edge fraction” denotes the leading-edge size divided by the gene-set size. 841 842 Knowledge graph creation and LLM inference 843 844 To identify unexplored relationships between our most explainable genes and endometriosis in a sparsely 845 characterised research area, we applied a literature -driven knowledge graph and large language model -846 based inference method to systematically integrate and reason over published evidence. This yielded a 847 literature-grounded representation of gene -endometriosis relationships suitable for potential hypothesis 848 generation. Each candidate gene was treated as an independent query entity and used to retrieve all PubMed-849 indexed articles in which it was mentioned, using open -source software DeepSearch 850 (https://ds4sd.github.io/). From the resulting corpus, we extracted paragraph-level text segments containing 851 co-occurrences of each entity and the endometriosis. These texts were processed using named entity 852 recognition and relation extraction to identify biomedical entities such as genes, diseases, proteins, 853 pathways, and used to identify the semantic relationships amongst them. All entities and relationships were 854 represented as entity-entity-relation triples, and integrated into a unified knowledge graph, with literature 855 provenance retained at document and paragraph levels. All inferred associations were constrained by the 856 underlying graph structure and traceable to one or more publications. To enable LLM-based graph-enabled 857 inference to identify higher -order or implicit relationships supported by the literature and interpret the 858 shared relational patterns in the compact KG, we supplied a compact JSON summary to a locally deployed 859 LLM (GPT-20b) executed using Ollama inference engine. We also calculated relation-style similarity and 860 shared-entity overlap for all our gene pairs (see below) and provided that to the LLM to aid inference. 861 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 862 Relation-style similarity (weighted Jaccard similarity) 863 864 To quantify similarity of connected entities and relations between the five important genes in the knowledge 865 graph, we calculated a weighted Jaccard similarity for each pair of important genes, quantifying the 866 similarity of their relation-type profiles across all shared entities. For each pair of important genes and for 867 each entity common to both genes, we calculated an intersection weight by summing the minimum 868 frequency observed for each relation type, as follows: 869 870 𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡𝑖𝑜𝑛 = J ! 𝑚𝑖𝑛K𝑐"(𝑟), 𝑐#(𝑟)L 871 872 Next, for that specific entity we captured the union weight which measures total combined relation usage, 873 so taking the maximum of the two counts as follows: 874 875 𝑢𝑛𝑖𝑜𝑛 = J ! 𝑚𝑎𝑥K𝑐"(𝑟), 𝑐#(𝑟)L 876 877 We then calculate the Jaccard similarity (intersection/union) with 1 indicating that the two genes share 878 exactly the same relation types (for shared entities) with same frequencies, and 0 indicating the two genes 879 share no common relation types for any shared entities. 880 881 Normalized shared entity overlap 882 883 For each gene pair, to account for substantial variation in the number of triples per gene due to varied 884 representation in the literature, we also calculated shared entity overlap directionally as the proportion of 885 shared entities from gene1 that were also present in gene2, and reciprocally from gene2 to gene 1. This 886 normalization allows overlap to be evaluated relative to each genes individual representation in the corpus. 887 888 Shared-entity overlap (Jaccard similarity) 889 890 To complement the relation -style similarity we also calculate a shared-entity Jaccard to quantify shared 891 entities between each gene pair. For each gene pair, this is defined as the proportion of shared (overlapping) 892 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint entities relative to the total number of distinct entities linked to each gene. For two genes 𝑔" and 𝑔#, the 893 shared-entity jaccard is calculated as follows: 894 895 𝐽$%&'&'$( (𝑔", 𝑔#) = |𝐸(𝑔") ∩ 𝐸(𝑔#)| |𝐸(𝑔") ∪ 𝐸(𝑔#)| 896 897 898 899 Declarations 900 901 Ethics approval and consent to participate 902 Not Applicable – all data used in this study was publicly available from the Gene Expression Omnibus 903 (GEO) 904 905 Consent for publication 906 Not Applicable 907 Availability of data and materials 908 All datasets used in this study are bulk RNA -seq data that is publicly available, Table 1 gives a 909 summary of the cohorts used in this study , detailing Gene Expression Omnibus (GEO) ID’s. The 910 foundation models and AutoXAI4Omics used in this study are all open source (see citations in 911 text). However, we describe a new interpretability method for the FM embeddings that we call 912 CA-IG. We will release this code as an open source repo (at https://github.com/IBM/) upon 913 acceptance of publication. 914 915 Competing Interests 916 The authors declare that they have no competing interests. Commercial affiliations do not alter 917 adherence of the authors to the journal policies on sharing data and materials. 918 919 Funding 920 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint This work was supported by the Hartree National Centre for Digital Innovation, a collaboration 921 between STFC and IBM funded by UKRI. 922 923 Authors’ Contributions 924 LJG & JK designed the study, LJG/JK/RT/MA/SC/AE contributed towards workflow conceptualization, 925 LJG, JK & NM performed the AI technical analysis, NM, JK & LJG performed the KG technical work with 926 contributions from JB, NC created the CA-IG methodology supervised by LJG & JK, JK, LJG & NM wrote 927 and edited the manuscript, MM contributed to BMFM embedding extraction, transcriptomic data collection 928 and curation, all authors reviewed code and/or manuscript. 929 930 931

932 This work was supported by the Hartree National Centre for Digital Innovation, a collaboration 933 between STFC and IBM. 934 935 936 937 938 939

940 941 1. Rogers PA, D’Hooghe TM, Fazleabas A, Gargett CE, Giudice LC, Montgomery GW, et al. 942 Priorities for endometriosis research: recommendations from an international consensus 943 workshop. Reproductive sciences. 2009;16(4):335–46. 944 2. Cox E. Dismissed, ignored and belittled”: the long road to endometriosis diagnosis in the 945 UK. Endometriosis UK. 2024; 946 3. Zondervan KT, Becker CM, Missmer SA. Endometriosis. New England Journal of Medicine. 947 2020;382(13). 948 4. Gete DG, Doust J, Mortlock S, Montgomery G, Mishra GD. Associations between 949 endometriosis and common symptoms: findings from the Australian Longitudinal Study on 950 Women’s Health. American Journal of Obstetrics and Gynecology. 2023 Nov;229(5):536.e1-951 536.e20. 952 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 5. Meuleman C, Vandenabeele B, Fieuws S, Spiessens C, Timmerman D, D’Hooghe T. High 953 prevalence of endometriosis in infertile women with normal ovulation and normospermic 954 partners. Fertility and sterility. 2009;92(1):68–74. 955 6. Simoens S, Dunselman G, Dirksen C, Hummelshoj L, Bokor A, Brandes I, et al. The burden 956 of endometriosis: costs and quality of life of women with endometriosis and treated in 957 referral centres. Human reproduction. 2012;27(5):1292–9. 958 7. Nisenblat V, Bossuyt PM, Farquhar C, Johnson N, Hull ML. Imaging modalities for the non-959 invasive diagnosis of endometriosis. Cochrane Database of Systematic Reviews. 2016;(2). 960 8. Nisenblat V, Bossuyt PM, Shaikh R, Farquhar C, Jordan V, Scheffers CS, et al. Blood 961 biomarkers for the non-invasive diagnosis of endometriosis. Cochrane Database of 962 Systematic Reviews. 1996;2016(5). 963 9. Moustafa S, Burn M, Mamillapalli R, Nematian S, Flores V, Taylor HS. Accurate diagnosis 964 of endometriosis using serum microRNAs. American journal of obstetrics and gynecology. 965 2020;223(4):557-e1. 966 10. Zhang H, Zhang H, Yang H, Shuid AN, Sandai D, Chen X. Machine learning-based 967 integrated identification of predictive combined diagnostic biomarkers for endometriosis. 968 Front Genet. 2023;14:1290036. 969 11. Kucukakcali Z, Akbulut S, Colak C. Prediction of genomic biomarkers for endometriosis 970 using the transcriptomic dataset. World Journal of Clinical Cases. 2025;13(20). 971 12. Akter S, Xu D, Nagel SC, Bromfield JJ, Pelch KE, Wilshire GB, et al. GenomeForest: An 972 Ensemble Machine Learning Classifier for Endometriosis. AMIA Jt Summits Transl Sci 973 Proc. 2020;2020:33–42. 974 13. Bendifallah S, Dabi Y, Suisse S, Jornea L, Bouteiller D, Touboul C, et al. MicroRNome 975 analysis generates a blood-based signature for endometriosis. Scientific Reports. 976 2022;12(1):4051. 977 14. Shi S, Huang C, Tang X, Liu H, Feng W, Chen C. Identification and verification of 978 diagnostic biomarkers for deep infiltrating endometriosis based on machine learning 979 algorithms. Journal of Biological Engineering. 2024;18(1):70. 980 15. Theodoris CV, Xiao L, Chopra A, Chaffin MD, Al Sayed ZR, Hill MC, et al. Transfer 981 learning enables predictions in network biology. Nature. 2023;618(7965):616–24. 982 16. Cui H, Wang C, Maan H, Pang K, Luo F, Duan N, et al. scGPT: toward building a 983 foundation model for single-cell multi-omics using generative AI. Nature methods. 984 2024;21(8):1470–80. 985 17. Hao M, Gong J, Zeng X, Liu C, Guo Y, Cheng X, et al. Large scale foundation model on 986 single-cell transcriptomics. bioRxiv. 2023; 987 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 18. Gélard M, Richard G, Pierrot T, Cournède PH. BulkRNABert: Cancer prognosis from bulk 988 RNA-seq based language models. bioRxiv. 2024;2024–06. 989 19. Dandala B, Danziger MM, Barkan E, Biswas T, Gurev V, Hu J, et al. BMFM-RNA: An 990 Open Framework for Building and Evaluating Transcriptomic Foundation Models. arXiv 991 preprint arXiv:250614861. 2025; 992 20. Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene expression and 993 hybridization array data repository. Nucleic acids research. 2002;30(1):207–10. 994 21. Lundberg SM, Lee SI. A unified approach to interpreting model predictions. Advances in 995 neural information processing systems. 2017;30. 996 22. Djos A, Martinsson T, Kogner P, Carén H. The RASSF gene family members RASSF5, 997 RASSF6 and RASSF7 show frequent DNA methylation in neuroblastoma. Mol Cancer. 2012 998 Dec;11(1):40. 999 23. Sangsin A, Kuptanon C, Srichomthong C, Pongpanich M, Suphapeetiporn K, Shotelersuk V. 1000 Two novel compound heterozygous BMP1 mutations in a patient with osteogenesis 1001 imperfecta: a case report. BMC Med Genet. 2017 Dec;18(1):25. 1002 24. Xu Z, Chen X, Song L, Yuan F, Yan Y. Matrix Remodeling-Associated Protein 8 as a Novel 1003 Indicator Contributing to Glioma Immune Response by Regulating Ferroptosis. Front 1004 Immunol. 2022 Feb 24;13:834595. 1005 25. Zhang M, Cai X, Liu J, Zhou J, Shi Q, Jiang Y, et al. A novel lncRNA lncSAMD11-1: 1 1006 interacts with PIP4K2A to promote endometrial decidualization by stabilizing FoxO1 1007 nuclear localization. The International Journal of Biochemistry & Cell Biology. 2022 1008 Oct;151:106280. 1009 26. Greenwood E, Maisel S, Ebertz D, Russ A, Pandey R, Schroeder J. Llgl1 prevents 1010 metaplastic survival driven by epidermal growth factor dependent migration. Oncotarget. 1011 2016 Sep 20;7(38):60776–92. 1012 27. Al-Hetty HRAK, Jabbar AD, Eremin VF, Jabbar AM, Jalil AT, Al-Dulimi AG, et al. The 1013 role of endoplasmic reticulum stress in endometriosis. Cell Stress and Chaperones. 2023 1014 Mar;28(2):145–50. 1015 28. Huang E, Zhang L, Lou J, Wang X, Chen L. Identified endoplasmic reticulum stress-related 1016 molecular cluster and immune characterization in endometriosis. Sci Rep. 2025 Nov 1017 4;15(1):38553. 1018 29. Choi J, Jo M, Lee E, Lee DY, Choi D. Dienogest regulates apoptosis, proliferation, and 1019 invasiveness of endometriotic cyst stromal cells via endoplasmic reticulum stress induction. 1020 Molecular Human Reproduction. 2020 Jan 1;26(1):30–9. 1021 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 30. Maytin EV, Ubeda M, Lin JC, Habener JF. Stress-Inducible Transcription Factor 1022 CHOP/gadd153 Induces Apoptosis in Mammalian Cells via p38 Kinase-Dependent and -1023 Independent Mechanisms. Experimental Cell Research. 2001 Jul;267(2):193–204. 1024 31. Ghosh AP, Klocke BJ, Ballestas ME, Roth KA. CHOP Potentially Co-Operates with 1025 FOXO3a in Neuronal Cells to Regulate PUMA and BIM Expression in Response to ER 1026 Stress. El-Rifai W, editor. PLoS ONE. 2012 Jun 28;7(6):e39586. 1027 32. Marola OJ, Syc-Mazurek SB, Libby RT. DDIT3 (CHOP) contributes to retinal ganglion cell 1028 somal loss but not axonal degeneration in DBA/2J mice. Cell Death Discov. 2019 Oct 1029 10;5(1):140. 1030 33. Yu X, Li W, Sun S, Li J. DDIT3 is associated with breast cancer prognosis and immune 1031 microenvironment: an integrative bioinformatic and immunohistochemical analysis. J 1032 Cancer. 2024;15(12):3873–89. 1033 34. Wang S, Hou P, Pan W, He W, He DC, Wang H, et al. DDIT3 Targets Innate Immunity via 1034 the DDIT3-OTUD1-MAVS Pathway To Promote Bovine Viral Diarrhea Virus Replication. 1035 Sandri-Goldin RM, editor. J Virol. 2021 Feb 24;95(6):e02351-20. 1036 35. Wang S, Hou P, Pan W, He W, He DC, Wang H, et al. DDIT3 Targets Innate Immunity via 1037 the DDIT3-OTUD1-MAVS Pathway To Promote Bovine Viral Diarrhea Virus Replication. 1038 Sandri-Goldin RM, editor. J Virol. 2021 Feb 24;95(6):e02351-20. 1039 36. Söderman J, Berglind L, Almer S. Gene Expression-Genotype Analysis Implicates GSDMA , 1040 GSDMB , and LRRC3C as Contributors to Inflammatory Bowel Disease Susceptibility. 1041 BioMed Research International. 2015;2015:1–10. 1042 37. Hodzic D, Kong C, Wainszelbaum MJ, Charron AJ, Su X, Stahl PD. TBC1D3, a hominoid 1043 oncoprotein, is encoded by a cluster of paralogues located on chromosome 17q12. 1044 Genomics. 2006 Dec;88(6):731–6. 1045 38. Wang B, Chen D, Hua H. TBC1D3 family is a prognostic biomarker and correlates with 1046 immune infiltration in kidney renal clear cell carcinoma. Molecular Therapy - Oncolytics. 1047 2021 Sep;22:528–38. 1048 39. Australian Cancer Study, Australian Ovarian Cancer Study Group, Pharoah PDP, Tsai YY, 1049 Ramus SJ, Phelan CM, et al. GWAS meta-analysis and replication identifies three new 1050 susceptibility loci for ovarian cancer. Nat Genet. 2013 Apr;45(4):362–70. 1051 40. Hodzic D, Kong C, Wainszelbaum MJ, Charron AJ, Su X, Stahl PD. TBC1D3, a hominoid 1052 oncoprotein, is encoded by a cluster of paralogues located on chromosome 17q12. 1053 Genomics. 2006 Dec;88(6):731–6. 1054 41. Guitart X, Porubsky D, Yoo D, Dougherty ML, Dishuck PC, Munson KM, et al. Independent 1055 expansion, selection and hypervariability of the TBC1D3 gene family in humans [Internet]. 1056 Genetics; 2024 [cited 2025 Dec 15]. Available from: 1057 http://biorxiv.org/lookup/doi/10.1101/2024.03.12.584650 1058 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 42. Shibel R, Sarfstein R, Nagaraj K, Lapkina-Gendler L, Laron Z, Dixit M, et al. The Olfactory 1059 Receptor Gene Product, OR5H2, Modulates Endometrial Cancer Cells Proliferation via 1060 Interaction with the IGF1 Signaling Pathway. Cells. 2021 Jun 12;10(6):1483. 1061 43. Qin S, Dorschner RA, Masini I, Lavoie‐Gagne O, Stahl PD, Costantini TW, et al. TBC1D3 1062 regulates the payload and biological activity of extracellular vesicles that mediate tissue 1063 repair. The FASEB Journal. 2019 May;33(5):6129–39. 1064 44. Wagner M, Hicks C, El-Omar E, Combes V, El-Assaad F. The Critical Role of Host and 1065 Bacterial Extracellular Vesicles in Endometriosis. Biomedicines. 2024 Nov 12;12(11):2585. 1066 45. Chen S, Birk DE. The regulatory roles of small leucine‐rich proteoglycans in extracellular 1067 matrix assembly. The FEBS Journal. 2013 May;280(10):2120–37. 1068 46. Woo J, Kwon SK, Kim E. The NGL family of leucine-rich repeat-containing synaptic 1069 adhesion molecules. Molecular and Cellular Neuroscience. 2009 Aug;42(1):1–10. 1070 47. Adachi M, Nasu K, Tsuno A, Yuge A, Kawano Y, Narahara H. Attachment to extracellular 1071 matrices is enhanced in human endometriotic stromal cells: a possible mechanism underlying 1072 the pathogenesis of endometriosis. European Journal of Obstetrics & Gynecology and 1073 Reproductive Biology. 2011 Mar;155(1):85–8. 1074 48. Mariadas H, Chen JH, Chen KH. The Molecular and Cellular Mechanisms of Endometriosis: 1075 From Basic Pathophysiology to Clinical Implications. IJMS. 2025 Mar 10;26(6):2458. 1076 49. Miller JE, Ahn SH, Marks RM, Monsanto SP, Fazleabas AT, Koti M, et al. IL-17A 1077 Modulates Peritoneal Macrophage Recruitment and M2 Polarization in Endometriosis. Front 1078 Immunol. 2020 Feb 14;11:108. 1079 50. Shi JL, Zheng ZM, Chen M, Shen HH, Li MQ, Shao J. IL-17: an important pathogenic factor 1080 in endometriosis. Int J Med Sci. 2022;19(4):769–78. 1081 51. Ahn SH, Edwards AK, Singh SS, Young SL, Lessey BA, Tayade C. IL-17A Contributes to 1082 the Pathogenesis of Endometriosis by Triggering Proinflammatory Cytokines and 1083 Angiogenic Growth Factors. The Journal of Immunology. 2015 Sep 15;195(6):2591–600. 1084 52. Zhang X, Xu H, Lin J, Qian Y, Deng L. Peritoneal fluid concentrations of interleukin‐17 1085 correlate with the severity of endometriosis and infertility of this disorder. BJOG. 2005 1086 Aug;112(8):1153–5. 1087 53. Delbandi AA, Mahmoudi M, Shervin A, Farhangnia P, Mohammadi T, Zarnani AH. 1088 Increased circulating T helper 17 (TH17) cells and endometrial tissue IL-17-producing cells 1089 in patients with endometriosis compared with non-endometriotic subjects. Reproductive 1090 Biology. 2025 Jun;25(2):101019. 1091 54. Hirata T, Osuga Y, Takamura M, Saito A, Hasegawa A, Koga K, et al. Interleukin-17F 1092 increases the secretion of interleukin-8 and the expression of cyclooxygenase 2 in 1093 endometriosis. Fertility and Sterility. 2011 Jul;96(1):113–7. 1094 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint 55. Podgaec S, Abrao MS, Dias JA, Rizzo LV, De Oliveira RM, Baracat EC. Endometriosis: an 1095 inflammatory disease with a Th2 immune response component. Human Reproduction. 2007 1096 May 1;22(5):1373–9. 1097 56. Wu MH, Yang BC, Hsu CC, Lee YC, Huang KE. The expression of soluble intercellular 1098 adhesion molecule-1 in endometriosis. Fertility and Sterility. 1998 Dec;70(6):1139–42. 1099 57. Keenan JA, Chen TT, Chadwell NL, Torry DS, Caudle MR. Interferon‐Gamma (IFN‐γ) and 1100 Interleukin‐6 (IL‐6) in Peritoneal Fluid and Macrophage‐Conditioned Media of Women With 1101 Endometriosis. American J Rep Immunol. 1994 Oct;32(3):180–3. 1102 58. Al-Hetty HRAK, Jabbar AD, Eremin VF, Jabbar AM, Jalil AT, Al-Dulimi AG, et al. The 1103 role of endoplasmic reticulum stress in endometriosis. Cell Stress and Chaperones. 2023 1104 Mar;28(2):145–50. 1105 59. Ekiz-Yilmaz T, Isildar B, Gezer A, Kankaya D, Cansiz-Ersoz C, Kayisli UA, et al. The role 1106 of unfolded protein response in the pathogenesis of endometriosis: contribution of peritoneal 1107 fluid. Reproductive BioMedicine Online. 2021 Jan;42(1):1–15. 1108 60. Tellez-Quijorna C, Juan-Lopez A, Eritja N, Llobet-Navas D, Devis-Jauregui L. Dysregulated 1109 autophagy in endometriosis: molecular mechanisms, controversies, and clinical implications. 1110 Reprod Biol Endocrinol. 2025 Dec 13;24(1):8. 1111 61. Mariadas H, Chen JH, Chen KH. The Molecular and Cellular Mechanisms of Endometriosis: 1112 From Basic Pathophysiology to Clinical Implications. IJMS. 2025 Mar 10;26(6):2458. 1113 62. Hao M, Gong J, Zeng X, Liu C, Guo Y, Cheng X, et al. Large scale foundation model on 1114 single-cell transcriptomics. bioRxiv. 2023; 1115 63. Strudwick J, Gardiner LJ, Denning-James K, Haiminen N, Evans A, Kelly J, et al. 1116 AutoXAI4Omics: an automated explainable AI tool for omics and tabular data. bioRxiv 1117 [Internet]. 2024; Available from: 1118 https://www.biorxiv.org/content/early/2024/03/28/2024.03.25.586460 1119 64. Freund Y, Schapire RE. A decision-theoretic generalization of on-line learning and an 1120 application to boosting. Journal of Computer and System Sciences. 1997;55(1):119–39. 1121 65. Sundararajan M, Taly A, Yan Q. Axiomatic attribution for deep networks. In: International 1122 conference on machine learning. PMLR; 2017. p. 3319–28. 1123 66. Fang Z, Liu X, Peltz G. GSEApy: a comprehensive package for performing gene set 1124 enrichment analysis in Python. Lu Z, editor. Bioinformatics. 2023 Jan 1;39(1):btac757. 1125 67. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a 1126 comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 1127 2016 Jul 8;44(W1):W90–7. 1128 1129 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted February 25, 2026. ; https://doi.org/10.64898/2026.02.24.707424doi: bioRxiv preprint