Unified Predictive Model for Endometriosis: Merging Clinical, Self-reporting and Genetic Information

UKBB data extraction and processing

The UK Biobank (UKBB) is a population-based database with detailed medical, genotyping, and lifestyle information covering 500k people at ages 40-69 at time of recruitment [29]. UKBB recruited the participants from 2006-2010 from across the UK. All analyses were based on the 2019 UKBB release.

We focused on Caucasian women by limiting the analysis to participants who self-reported themselves as British, Irish, or other “white” background [codes 1, 1001, 1002, 1003, respectively, in Ethnic background, UKB data field 21000]) and classified as Caucasians based on their genetic ancestry (Genetic ethnic group, data-field 22006). We further focused on participants evaluated at age 40-70 (dated in 2006) and removed genetic relatives, by keeping only one representative of each kinship group of related individuals. This resulted in a dataset that includes 145,671 participants. Disease classification is based on clinical information encoded by ICD-10 codes. We used the main or secondary diagnosis (UKBB data fields 41202 and 41204, respectively) with the age of the diagnosis. We addressed each data field according to the missing information included. The fraction of missing data is restricted to missing measurements of participants that did not know the answer (e.g., breastfed as a baby). Notably, for some cases, the information is only relevant for a subset of the studied population. In cases where multiple values were reported for a specific field (e.g., BMI from repeated visits), only the last value was considered. Data fields that were found related to endometriosis by the literature (and consulting with clinicians) were collected, along with all of the participants’ documented ICD-10 code diagnosis.

A protocol for age-dependent matching of endo group and control group was by performing a stochastic matching process between the two groups. The objective of this protocol is to keep the majority of the samples while matching the year of birth distribution. In practice, we randomly choose 71,088 samples from the group of women without endometriosis diagnosis (control group) in an amount that imitates the year of birth distribution of women with endometriosis diagnosis (endo group). The rest of the analysis was performed on the yearly matched set. See Supplementary Text S1 for the pseudocode used.

Genetic Analysis

The UKBB released genotyped data for all participants. The genotyping scheme is based on 805,426 preselected genetic variations. Based on the imputation protocol, the number of variants is expanded to about 9 M variants passed quality control [30]. We used the OpenTargets (OT) platform to select current knowledge on endometriosis genetics [31]. OT is a public database that unifies evidence for drugs, their targets, and their associations with human diseases. We used the genetic platform that compiled the top-scored variants from GWAS summary statistics as extracted from the GWAS catalog [32]. We used the OT genetic association scoring system to extract an informative list of variants associated with endometriosis. We gathered a list of 189 SNPs from all 221 genes from OT (based on OT quality criteria, some genes lack associated SNPs). We extracted the SNPs associated with endometriosis as reported by OT. A total of 65 unique genetic variants were used in our model.

Machine learning methodology

We tested several models including Random forest, Logistic regression, Linear discriminant analysis, and compared their performance. We also applied CatBoost that belong to a family of trees-based gradient boosting algorithms that perform well in big data with missing data [33,34]. In a nutshell, in each step of the algorithm, a decision tree base learner is created, using the previous iterations’ decision tree residuals as a gradient for minimizing the current tree’s loss-function. For each iteration, CatBoost uses a random permutation of the training set. The subset is used in order to build the decision tree and to build target statistics for the categorical features by mapping these features into a continuous space [35]. We trained three types of models according to the type of data used for training: (a) Attributes and measurements that were compiled from the reported risk factor for endometriosis in literature and other fields that were proposed by medical experts (Supplementary Table S1). (b) Medical diagnoses, as indexed by ICD-10 codes. (c) Genetic variants based on endometriosis GWAS from UKBB marker SNPs (Supplementary Table S2).

We used SHAP (SHapley Additive exPlanations) to estimate the features’ importance [34]. SHAP values give a numerical estimate of the marginal impact of a feature, given all other features.

Feature engineering

In addition to the UKBB data fields, we engineered features which were not explicitly found in the UKBB. Estrogen exposure, for example, was calculated by reducing the age of menarche from the age of menopause. Many of the features from the ICD-10 diagnosis fields were extracted from the UKBB and converted prior to their use in the predictive model (Supplementary Table S3). We calculated from the reported dates of any diagnosis available in the UKBB the age when the participant was diagnosed for each of the ICD-10 available for that person. A feature of the amount of ICD-10 diagnoses was calculated by summing the diagnoses available in the medical record that were accumulated prior to endometriosis diagnosis age. In this case, for the control group, a matching protocol was performed in order to determine the age threshold for such counting.

Statistical tests

We applied post-hoc univariable analysis using Kruskal-Wallis for continuous variables and Pearson’s chi-squared test for binary variables. For each feature, we calculated the standardized mean difference (SMD) as its summary statistics. SMD expresses the size of the effect relative to the variability observed. Formally we measured the mean outcome between endometriosis patients and the control group relative to the standard deviation of the outcome among control participants. The univariable analysis was limited to Q1-Q3 to improve statistical robustness.

Unification of data from UKBB: Case-control population-based groups

The primary goal of this study was to review current risk factor knowledge and evaluate its contribution to endometriosis prediction. To this end, we systematically collected a set of phenotypes and measurements extracted from the UKBB database. As a population-based resource, the UKBB is based on standardized data collection protocols. The UKBB includes over 500,000 participants collected from 23 medical centers across the UK, who were recruited over the years 2006–2010 for participants aged 49–70. We have retrospectively analyzed personalized clinical information on diagnosis, medical procedures, lifestyle, personal genetics, self-reporting, and nurse interview reports. Following strict filtration steps (see Methods), we analyzed 148,571 women, among whom 5924 were diagnosed with endometriosis (ICD-10: N80, Data field).

Table 1 lists a selected sample for the different data types (e.g., physical measurement) that were used in this study. Note that the extracted UKBB fields cover information that is binary, contentious, or divided to discrete categories. The data extraction following the filtration scheme covered 970 diagnoses, 65 genetic variants, 46 life-style and physical measures. Supplementary Table S1 lists all the life-style and physical measures UKBB data fields extracted and the degree covered by the 148.5k women included in this study. The extraction of data was motivated by endometriosis risk factors previously studied and expanded according to an input from medical experts.

The data of the UKBB was obtained by the participants’ medical records or by questioners and exams at assessment centers. Despite the effort to standardize and fill all data fields, listed in Supplementary Table S1, some attributes and measurements suffer from a substantial fraction of missingness. For example, only 2.7% of cases lack menarche age, while the ages of the first and last age of depression episodes are missing for 78.5% of cases. Figure 2 lists the variables (Supplemental Table S1) for covering missing data at a range of 2% to 80%. Note that the fraction of missing data is calculated from the number of participants that were diagnosed with the relevant diagnosis. For example, the ‘age of the first episode of depression’ is only valid to those who replied positively to ‘ever felt depression’. Among those subjects, 80% had not reported on the age at the first episode of depression.

• Download figure
• Open in new tab

Figure 1.

Ranked list of variables with the percentage of the missing data. A full list of extracted attributes is available in Supplementary Table S1.

• Download figure
• Open in new tab

Figure 2.

Univariate analysis for endometriosis. A ranked list of attributes (total 44) associated with endometriosis diagnosed and control groups by the standardized mean difference (SMD). SMD values <−0.2 and >0.2 are colored orange to indicate the attribute with a substantial effect size. The statistics were based on the median calculated for the Q1-Q3 values. The asterisk next to the description of the attribute is the case with p-value <0.05 for univariate tests of case and control (see Methods). For a univariate statistical test and results, see Supplementary Table S1.

Univariate statistics of control and endometriosis patients from the UKBB

A post hoc statistical test was performed to assess the contribution of each individual measurement. Numerous attributes have been previously reported as risk factors for endometriosis. Figure 2 shows the differences between the endo group and the control group based on SMD (see Methods). Each attribute was independently analyzed by including the median values (Q1, Q3) and calculating the statistical significance of its effect size (Supplementary Table S1). Setting the SMD threshold at 0.2, only 6 (out of 44) attributes are strongly associated with risk for endometriosis. The number of live births and the age at cancer and diabetes diagnosis (UKBB fields of 2734 and 40008, respectively) suggest a lower risk for endometriosis. The most significant variable in accordance with an increased risk of endometriosis is the year of birth (SMD 0.44) followed by irritable bowel syndrome (IBS). The rest of the measurements had smaller effect sizes. For detailed information, see Supplemental Table S1.

The calculated effect sizes associated with most of the attributes associated with endometriosis (e.g., menarche age, BMI, height, birth weight) were low. Other attributes failed to meet statistical significance (e.g., smoking, height, coffee consumed). The list shown in Figure 2 did not overlap with known risk factors for endometriosis as reported in the literature.

Endometriosis is a complex condition and assessing the risk according to the assessment of each attribute independently of the others cannot capture the interactions and the non-additive contributions of specific factors. A likely scenario is that different factors (each carrying a marginal effect) interact, and their combination provides valuable predicting power. Moreover, the extracted and engineered features belong to multiple types. Some attributes are continuous (e.g., BMI), others are binary (e.g., having a specific ICD-10) and many are assigned categories (e.g., smoking habits). Thus, we seek a method that considers any variable irrespectively of its type. For the goal of developing a predictive model for endometriosis, we applied a multivariate machine learning-based framework. A scheme of the analyses and processes for creating a predictive model for endometriosis using the UKBB data is shown (Figure 3). In brief, following filtration (see Methods), a screening process was applied (see Methods), resulting with 148,571 participants, out of whom 5,924 were diagnosed with endometriosis. The data were split to disjoint 80% training and 20% test sets. We further analyze the data and its distribution to account for internal year-dependent biases (Figure 3, Data processing).

Figure 4 shows the distribution of the participants in the study for women that were not diagnosed (control group) and those diagnosed with endometriosis (endo group). There was a significant difference in the year of birth distribution among women with and without endometriosis (U-test p-value 2.2e-239). Evidently, with very significant statistical differences, it is anticipated that a bias by the year of birth for the endo group is probably a reflection of establishing the diagnosis protocol and increasing awareness. To overcome this bias, we created a matched set for each year to cancel out the original year of birth differences. Repeating the U-test after applying the matching protocol resulted in insignificant difference between the control group and the endo group. The rest of the analysis was performed on the age-matched data.

• Download figure
• Open in new tab

Figure 4.

The distribution of the control and endo-groups along the year of birth (Left). Following a protocol for yearly matching schemes, the bias was removed. And each year a matched proportion of control and endo-groups remains stable throughout (for detailed protocol, see Supplementary Text S1).

Predictive risk model for endometriosis

We used the receiver operating characteristic area under the curve (roc-AUC) as the evaluation metric. The CatBoost model trained for 1000 iterations using early stopping on a separate held out validation subset. After a screening process (Figure 3), the data was separated into three main categories according to the type of data used for training. These categories are the basis for three models labelled a, b and c according to the type of data used as input (see Methods): (a) Attributes and measurements form UKBB (Figure 3, Supplementary Table S1); (b) Medical diagnoses, as indexed by ICD-10 codes (Supplementary Tables S3); and (c) Genetic variants based on endometriosis GWAS from UKBB marker SNPs (Supplementary Tables S2).

In preparation for model b, we first tested whether differences between the cases and controls could be derived from the associated vector of ICD-10 diagnoses (UKBB data fields 130000–132606). These UKBB data fields provide the dates of the participants’ initial appearance of any reported medical diagnosis. The dates were converted into the age of diagnosis for each woman. We asked whether the set of diagnoses is informative for endometriosis prediction. The rationale is to assess whether other diagnoses preceding the definitive endometriosis diagnosis, carry a predictive power towards endometriosis. For each participant in the control group, a threshold age for the diagnosis masking was randomly chosen from the endometriosis diagnosis age, such that the threshold distribution in the control group is equal to the distribution of endometriosis diagnosis age. The median number of diagnoses prior to that of endometriosis for the control and endo-group was 1 and 4, respectively (Figure 5A).

• Download figure
• Open in new tab

Figure 5.

ICD-10 in control and endo-groups. (A) The distribution of the amount of ICD-10 diagnoses in the control and endo-groups (orange and blue, respectively) was significant using Mann-Whitney U-test (p-value <0.001) and SMD = 0.471. The median value of the number of ICD-10 diagnoses per individual for the control and endo-groups is 1 and 4 respectively. (B) Partition of all 222 all statistically significant informative features from the ICD-10 based model (U-test, p value <0.05). Each feature was tested for the statistical difference of the control and the endo-group. The partition is according to the ICD-10 level1 first letter (A-Q). The level 1 letters with less than 10 features are unified (‘others’). (C) Ranked list of the top 40 ICD-10 that statistically differentiate ranked by the p-value <1e-11. These 40 ICD-10 are color coded as in B by level 1 ICD-10 index. The detailed information of the feature and its ICD-10 level 4 information is available in Supplementary Table S3.

Supplementary Table S3 shows the percentage of ICD-10 terms associated with women with and without endometriosis for 755 age associated diagnoses (see Methods). While only 7% of the control group have >10 ICD-10 diagnoses, there are 11% of the endo group with more than 30 ICD-10 diagnoses. Each age-converted ICD-10 was tested for the statistical difference of the control and the endo-group. For 222 items the “age of first reported diagnosis” resulted in p-value <0.05 by non-parametric statistical test (Supplementary Table S3). Figure 5B shows the partition of these 222 items according to ICD-10 indexing method (level-1; marked as A to Q; Supplementary Table S4). The abundant ICD-10 level 1 includes diseases of the genitourinary system (N) followed by diseases of the digestive system (K), diseases of the musculoskeletal system and connective tissue (M) and diseases of the respiratory system (J). The significant of diseases of the respiratory system (J) and viral and parasite infection (B) is less evident.

Figure 5C shows a ranked list of the most significant ICD-10 items according to U-test statistical results with pelvic and genital organs that prevail. Specifically, the most significant ICD-10 items include N73 (pelvic inflammatory diseases), N81 (female genital prolapse), noninflammatory disorders of ovary, fallopian tube, and broad ligament (N83) and of uterus, except cervix (N85), polyp of female genital tract (N84) and excessive, frequent and irregular menstruation (N92). While the knowledge on endometriosis corroborate the relevance of diseases associated with N, K and M, and to a lesser level of significance also diseases of the respiratory system (J).

In preparation for the machine learning predictive genetic model (model c), we collected variants from GWAS of endometriosis as an input for the predictive model. A list of 65 genetic variants associated with 35 different genes was compiled from11 major publications including large meta analyses [26] (Detailed in Supplementary Table S2).

Figure 6 shows the results from the AUC and the ROC curve for 5 models that are based on the major data type categories (marked a, b, c; see Methods) and their combination. The predictive models for each of the data types (a-c) and their combinations are shown for the combination of recall and precision (Figure 6A) and the rec-AUC of all 5 models in Figure 6B. Developing a model based on the 65 variants from the GWAS catalog (model c) indicated that training the model on genotypic data resulted in roc-AUC of 0.52. A recent population-based polygenic risk score (PRS) analysis for endometriosis showed only 2-3% of the variance explained by the SNPs [36], consistent with the modest improvement in the performance of model c. We found that model c, in combination with model a (measurements and attributes from self-reporting and lifestyle data) and model b (age-converted ICD-10 for diagnoses prior to endometriosis) resulted in higher AUC (increased from 0.77 to 0.81; Supplementary Table S5).

• Download figure
• Open in new tab

Figure 6.

Performances of the prediction models for endometriosis. (A) Precision-Recall curves for 5 CatBoost models. The models differ by training data with the UKBB attributes and measurements (model a), the collection of the ICD-10 prior to endometriosis diagnosis age (model b), and the GWAS of endometriosis genetic variants (model c). A combination of training data of a & b and a combined model that includes a, b & c. (B) ROC curves for the same set of 5 models as in A. The diagonal line marks a random no discrimination line (AUC = 0.5). (C) Comparison of the roc-AUC of 5 different algorithms for the combined set of features input a, b & c. XGBoost and CatBoost resulted in the highest roc-AUCs.

We repeated training with input a, b and c to test the performance of additional machine learning models (Figure 6C). The results of the models performed by Random forest, Logistic regression, Linear discriminant analysis, XGBoost and CatBoost algorithms are shown. The CatBoost algorithm of the combined model outperformed other models, followed by XGBoost (Supplementary Table S6). A much smaller AUC was associated with algorithms including K-nearest neighbors (KNN), Naive Bayes (NB) and support vector machines (SVM).

Informative features and interpretability of the combined model

We further evaluated the contribution of each feature on the combined model that was trained on 3 groups of features using SHAP, an explainable AI tool. Figure 7 shows the top 20 features ranked by SHAP. About a third of these features are associated with features of the age-dependent ICD-10, level 1 (Figure 5B), with the rest derived from the features associated with measurements and UKBB attributes. The top features are the length of the menstrual cycle and the age of the first live birth. Note that none of the genetic variants was selected to be informative among the top 20 features. Figure 7 also emphasizes the limited overlap between SHAP informative features and the attributes with significant as significant SMD from the univariate test (Figure 2).

• Download figure
• Open in new tab

Figure 7.

Top 20 features from the combined model using SHAP (an explainable AI tool). Variables are ranked in descending order of their SHAP value. The values reported show the contribution of each of the features according to the impact of that feature on the model outcome (i.e., endometriosis). Each dot in the plot represents a subject patient’s feature value for that variable (vertical axis). Color reflects the scale of the feature’s value. Color shows whether that variable is high (red) or low (blue) for that observation, gray depicts no data or a categorical feature.

The significant SHAP values supports the contribution of noninflammatory disorders of ovary, fallopian tube and broad ligament (SHAP value of 0.134), and excessive, frequent and irregular menstruation (N-92, SHAP value of 0.124). The informative features ranked by SHAP (e.g., estrogen exposure, reports of IBS) also displayed a strong deviation in occurrence in the endo group and control groups (Supplementary Table S3). However, statistically significant features from the ICD-10 diagnoses by age abundant in the endo group relative to control group not selected as informative features by SHAP. This list includes the age of the first occurrence in N39 (other disorders of urinary system), I10 (essential, primary, hypertension) and D50 (iron deficiency anemia) with p-values of 7e-55, 6e-42 and 2e-35, respectively.

Model’s limitation

Almost all analyzed data used for our models are based on measurements observed from women after their menopause age. Thus, the up-to-date diagnostic measurements are unavailable. The presented models (Figure 6C) are not designed as tools for diagnosis. However, we engineered features that include information collected prior to the date of diagnosis of endometriosis. As diagnosis is likely to be delayed, the partition before and after endometriosis diagnosis might be inaccurate. Another aspect that may limit the generality of our model concern an unavoidable enrichment in women with symptomatic, severe endometriosis. We anticipate that data analyzed from these women may not represent mild manifestations of endometriosis. In terms of UKBB data quality, data fields of UKBB diagnosis that lack a time stamp cannot be determined whether occurred before endometriosis diagnosis.