Article Text
Abstract
Objectives To validate and test the generalisability of the SASKit-ML pipeline, a prepublished feature selection and machine learning pipeline for the prediction of health deterioration after a stroke or pancreatic adenocarcinoma event, by using it to identify biomarkers of health deterioration in chronic disease.
Design This is a validation study using a predefined protocol applied to multiple publicly available datasets, including longitudinal data from cohorts with type 2 diabetes (T2D), inflammatory bowel disease (IBD), rheumatoid arthritis (RA) and various cancers. The datasets were chosen to mimic as closely as possible the SASKit cohort, a prospective, longitudinal cohort study.
Data sources Public data were used from the T2D (77 patients with potential pre-diabetes and 18 controls) and IBD (49 patients with IBD and 12 controls) branches of the Human Microbiome Project (HMP), RA Map (RA-MAP, 92 patients with RA, 22 controls) and The Cancer Genome Atlas (TCGA, 16 cancers).
Methods Data integration steps were performed in accordance with the prepublished study protocol, generating features to predict disease outcomes using 10-fold cross-validated random survival forests.
Outcome measures Health deterioration was assessed using disease-specific clinical markers and endpoints across different cohorts. In the HMP-T2D cohort, the worsening of glycated haemoglobin (HbA1c) levels (5.7% or more HbA1c in the blood), fasting plasma glucose (at least 100 mg/dL) and oral glucose tolerance test (at least 140) results were considered. For the HMP-IBD cohort, a worsening by at least 3 points of a disease-specific severity measure, the "Simple Clinical Colitis Activity Index" or "Harvey-Bradshaw Index" indicated an event. For the RA-MAP cohort, the outcome was defined as the worsening of the "Disease Activity Score 28" or "Simple Disease Activity Index" by at least five points, or the worsening of the "Health Assessment Questionnaire" score or an increase in the number of swollen/tender joints were evaluated. Finally, the outcome for all TCGA datasets was the progression-free interval.
Results Models for the prediction of health deterioration in T2D, IBD, RA and 16 cancers were produced. The T2D (C-index of 0.633 and Integrated Brier Score (IBS) of 0.107) and the RA (C-index of 0.654 and IBS of 0.150) models were modestly predictive. The IBD model was uninformative. TCGA models tended towards modest predictive power.
Conclusions The SASKit-ML pipeline produces informative and useful features with the power to predict health deterioration in a variety of diseases and cancers; however, this performance is disease-dependent.
- DIABETES & ENDOCRINOLOGY
- Health informatics
- Inflammatory bowel disease
- Rheumatology
Data availability statement
Data are available in a public, open access repository. Data are available on reasonable request. The data used in this study come from the publicly available datasets of the HMP, RA-MAP and TCGA. R scripts and preprocessed data to replicate the study are provided at https://github.com/danrpalmer/saskit-ml.
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See: http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
STRENGTHS AND LIMITATIONS OF THIS STUDY
Prepublished protocol reduces the chance of false positive results.
Chosen feature integration methods are inherently explainable, avoiding ‘black box’ machine learning.
Comparatively low sample sizes in many of the analysed cohorts.
Following a prepublished protocol restricts the use of newer methods.
Introduction
The SASKit project is an observational study seeking to identify and evaluate biomarker signatures for the prediction of health deterioration, and ultimately survival, in patients that have suffered a (pancreatic) adenocarcinoma or stroke event, with a view to optimising the provision of care to these patients. In brief, ~50 stroke patients (already recruited at the time of writing), ~50 (pancreatic) adenocarcinoma patients and ~50 age-matched controls (drawn primarily from the spouses of patients in the study) are being recruited and observed over the duration of the project. Baseline measurements of numerous candidate biomarkers (gene and protein expression and a selection of clinical variables) are taken at recruitment, with clinical biomarkers also being measured at multiple further time points, and shall then be analysed and used to establish effective biomarker signatures of health deterioration.
We published the proposed analysis protocol for the SASKit endpoint analyses in 2020.1 The protocol begins with two independently executed feature integration steps: Gene Set Variation Analysis (GSVA)2 and ExprEssence.3 Subsequently, biostatistical analyses (based on Cox regression analysis4 and machine learning analyses (based on random survival forests5) are performed to build predictive models of health deterioration. Here, we evaluate the random survival forest analysis using, insofar as is possible, the features proposed in the SASKit protocol.1 The main deviations in this analysis are in the diseases studied (as one goal of this analysis is to test the generalisability of the protocol to other diseases), the sample sizes which for some datasets (notably Human Microbiome Project (HMP)-inflammatory bowel disease (IBD)) are very similar to the SASKit cohort but which for others (in particular the larger The Cancer Genome Atlas (TCGA) datasets) are very different, and the specific tissue types sampled (again, the TCGA datasets being the largest deviation). No deviations in the modelling strategy were made for this analysis, with the exception that where a particular feature specified in the SASKit protocol was not available or able to be calculated, it was not used.
We fixed the data analysis protocol before any data were generated;1 with any analyses based on other methods labelled as ‘exploratory’. Naturally, new methods are published on an ongoing basis, sometimes shown to outcompete the methods we selected in 2020. Numerous feature selection methods for time-to-event analysis are available and choosing the correct approach is a challenge,6 especially if a heterogeneous set of data is provided as input, for example, omics and non-omics data. The choice is dependent on the variety and amount of input data, as well as the target output. Further, the application domain is of relevance, as is the experimental design behind the study that delivers the data to be analysed. For instance, long-term follow up of cancer patients could be expected to feature stronger signals than short-term studies of ageing in healthy people. With regard to the methods used here, a 2021 (available online from March 2020) benchmarking paper recommended GSVA, particularly for situations where single sample scores were required,7 while in 2020, it was shown that random survival forest methods can handle multiomics data better than several other survival analysis methods.8 No recent benchmarking study has included ExprEssence.
Additionally, SASKit was designed to investigate the hypothesis that the health deterioration associated with disease may be mediated by cellular senescence, and thus, that cellular senescence-related processes will be detectable in the features contributing to the best predictor(s). It is increasingly apparent that senescent cells are linked as contributing or causal factors in several diseases.9 10 This association leads to the possibility of using the measurement of senescence markers or associated pathways as biomarkers for a variety of clinical outcomes. Circulating senescence-associated secretory phenotype (SASP) proteins, for example, have been associated with frailty and adverse postsurgery outcomes11 while the modulation of SASP gene expression in skin biopsies has been put forward as a predictor of clinical response in systemic sclerosis.12
Here, we evaluate the SASKit-ML pipeline, as defined in 20201 in publicly available datasets with sufficient similarity to the datasets used for the SASKit study design. This approach gives an idea of the generalisability of our method and suggests future refinements. Here, SASKit-ML refers specifically to the machine learning pipeline in the published protocol, as opposed to the biostatistics pipeline outlined in the same paper. By trying to match the data types and protocol as closely as possible to the canonical SASKit study, we aim to both assess the validity of the SASKit-ML pipeline (under the assumption that good performance in similar datasets would predicate good performance in the SASKit cohort) as well as its generalisability to diseases other than stroke and pancreatic ductal adenocarcinoma (PDAC).
Methods
Datasets
We found the NIH HMP13 14 dataset to be one of the few publicly available longitudinal datasets that includes high-throughput transcriptomics and proteomics, a wide range of clinical blood markers, as well as time-to-event endpoint data. The two branches of the HMP most appropriate for analysis by the SASKit-ML pipeline (detailed below) are the type 2 diabetes (T2D) and IBD cohorts, however, each presents its own challenges. The T2D cohort of the HMP is the most comparable, containing peripheral blood mononuclear cells (PBMC) longitudinal transcriptomics and plasma proteomics (whereas in the IBD cohort, the proteomics are only from stool samples), as well as various clinical data and markers and the clinically defined endpoints of onset of pre-diabetes or T2D itself. Both cohorts have controls based on patients who were recruited to the project based on the relevant suitability metrics (ie, at risk of pre-diabetes or presenting with symptoms of IBD) but who then did not meet the criteria under closer observation. While these controls can be used for the purposes of SASKit feature integration, they do not follow the patient versus healthy control dichotomy as clearly as those in the SASKit study, and so the signal could be expected to be lower in these datasets. Then again, there is strong literature evidence to support the role of cellular senescence in both T2D15 and IBD.16
Another suitable dataset is the Rheumatoid Arthritis Map (RA-MAP),17 a large database of RA molecular data that include an extensive, longitudinal study including transcriptomics, proteomics and clinical values for patients with RA.18 This study very closely matches SASKit, having both whole blood transcriptomics and plasma proteomics at baseline (compared with PBMC transcriptomics and serum proteomics in SASKit). Finally, RA-MAP also includes suitable controls, in the form of baseline data from a simultaneous vaccine study (ie, omics data measured in healthy people before the intervention was administered). Despite this dataset containing high-sensitive C reactive protein (hsCRP) data, as in the SASKit cohort, hsCRP was not considered as a feature due to its use in the definition of the endpoint (hsCRP is a component of "Simple Disease Activity Index (SDAI)).
Additionally, TCGA data19 can be viewed in relation to the (pancreatic) adenocarcinoma branch of the SASKit cohort, though with some important differences. TCGA data are composed of high-throughput transcriptomics from tumour tissue rather than blood. While TCGA does contain proteomics, the lack of suitable and sufficient control proteomics samples limits their use in the SASKit-ML pipeline. A final difference is that the control tissue samples from TCGA are taken from healthy tissue, adjacent to the tumour tissue and are thus paired. We used TCGA data purely to test the effectiveness of the SASKit transcriptomics feature integration approach. In total, 16 TCGA datasets were found to be appropriate for the SASKit-ML pipeline. The datasets used are outlined in online supplemental table S1.
Supplemental material
Aside from the diseases being studied, and the specific clinical markers and omics available for each dataset, a final and important difference between all of these datasets and the SASKit cohort is that the SASKit examination times and thus potential endpoint times are relatively limited (with exams at 0, 3, 12, 24 and 36 months). The T2D dataset on the other hand has 41 unique event times, the IBD dataset has 23 and the RA-MAP, most similar to SASKit in this regard, has 4. TCGA likewise has a typically high and varied number of endpoint times compared with the SASKit cohort. In all cases, missing data were handled by median imputation. Further details pertaining to the key elements of each cohort, including study dates, follow-up and eligibility criteria can be found in their respective publications.
Biomarker discovery (SASKit-ML pipeline)
The SASKit data analysis plan has been published in advance.1 The underlying hypothesis is that, based on clinical markers, gene expression-based pathway scores, subnetwork activity scores and protein expression measurements, clinically relevant predictors of health deterioration/survival can be learnt for (pancreatic) adenocarcinoma and ischaemic stroke. The biggest challenge in adapting our proposed methodology to other datasets is that sufficiently similar datasets are currently rare. As such, for this study, it was considered sufficient for a dataset to have, at least, RNAseq transcriptomics and time-to-event data available so that the transcriptomics feature integration steps could be tested. A summary of the data types used for each dataset can be seen in online supplemental table S1 and code for the reproduction of this analysis along with a detailed breakdown of the various steps can be found in the project GitHub repository (https://github.com/danrpalmer/saskit-ml).
Dataset summary and endpoints
To align each cohort to the SASKit methodology, we defined suitable endpoints for them along similar lines to those used in SASKit, that is, endpoints that represented a worsening of the studied disease and/or decline in quality of life. For the T2D study (59 patients at risk of pre-diabetes and 13 healthy controls), the endpoint was set as a worsening of the glycated haemoglobin (HbA1c) value to prediabetic levels (5.7% or more HbA1c in whole blood), a fasting plasma glucose of at least 100 mg/dL and/or an oral glucose tolerance test result of at least 140 mg/dL, in accordance with the thresholds of the American Diabetes Association.20 For the IBD study (consisting of 66 patients with IBD and 20 non-IBD cases), the endpoint was set as a worsening by at least 3 points of the disease-specific severity measure ("Simple Clinical Colitis Activity Index" (SCCAI) for ulcerative colitis patients and "Harvey-Bradshaw Index" (HBI) for Crohn’s disease patients) when compared with the baseline measurement (in line with Allegretti, et al.).21 Although SCCAI and HBI are defined for two different diseases, we expected that the blood-based markers for both diseases may overlap. Indeed, hsCRP is a general marker of inflammation and as such can be predictive of disease activity in both diseases.22 As such, ulcerative colitis and Crohn’s disease patient data were combined for the analysis. For the RA-MAP study (99 patients with RA, 22 healthy controls), the endpoint was reached if any of the following commonly used indicators23 of RA disease severity became true relative to the baseline: a worsening of the "Disease Activity Score 28" (DAS28) score or SDAI score by at least 5 points, and/or a worsening of the Health Assessment Questionnaire score by at least 1 point and/or an increase in the number of swollen and/or tender joints by at least 1. For the TCGA data, the endpoint was disease progression indicated by the progression-free interval as defined in Liu et al.24
Dataset cleaning and preprocessing
The SASKit-ML pipeline is designed to take raw count (transcriptomics) and raw expression values (depending on the proteomics technology used). As such it was important to align the data available for the validation datasets as closely as possible. For all datasets other than HMP-T2D (normalised by DESeq2) and RA-MAP (expression array), raw count transcriptomics were available. For the HMP-T2D dataset, the required samples were identified by taking as baseline the earliest time point with all required data types measured, relative to the study starting date of 3 January 2013. This does mean that for some patients, what we considered the baseline measurement for this analysis is not close to their entry point for the study (ie, in cases where there were numerous visits where only clinical data were collected before any samples were subjected to omics analysis). However, for most patients, there remained a significant amount of follow-up from this point, with a median time of 643 days from the baseline measurement to the final visit for each patient.
Feature generation
As specified in the prepublished data analysis plan,1 a maximum of 25 integrated features was considered for each dataset. Additionally, as per the data analysis plan, age and sex were included as covariates in the calculation of the strongest GSVA features, in addition to being included in the models themselves. Thus, the following features were generated for each dataset:
Up to five curated clinical markers.
Up to five transcriptomic GSVA features.
Up to five transcriptomic ExprEssence features.
If high throughput proteomics is available, these were supplemented by:
Up to five proteomic GSVA features.
Up to five proteomic ExprEssence features.
The details of how each of these features were chosen and generated are as follows.
Clinical markers
Only partial overlap with the SASKit clinical markers was possible for each dataset, and while these were chosen for their relevance to pancreatic adenocarcinoma and ischaemic stroke, they are all also general ageing/senescence-associated markers and so should also hold predictive power for other diseases. Specifically, the HMP-T2D dataset included neutrophil–lymphocyte ratio, hsCRP and albumin, the HMP-IBD dataset included only hsCRP, while the RA-MAP dataset included neutrophil–lymphocyte ratio.
Gene Set Variation Analysis
GSVA is a popular method of data integration that condenses expression level data into pathway activation scores, which can then be compared between conditions.2 In short, GSVA calculates an expression-level statistic for each gene in each sample using a kernel estimation of the cumulative density function, then uses these statistics to rank the genes within each sample. On this basis, a Kolmogorov-Smirnov-like rank statistic is calculated for each gene set, and a pathway activation score is given by calculating the difference between the largest and smallest values of the random walk, analogous to the Kuiper test statistic.
Here, for both transcriptomics and proteomics (where available), we calculated the GSVA pathway activation scores for the 182 KEGG pathways taken from the curated Human MSigDB C2 gene set.25 26 These KEGG pathway activation scores were then subject to a differential activation analysis, comparing the case to control samples for each dataset. As per the study protocol, the differential activation analysis was carried out using the limma R package, employing a linear model with age and sex as additional covariates. Pathways were then ranked by significance, or, in the case that no significant pathways were detected, ranked by fold-change instead and the top-5 most strongly changing pathways taken as features.
ExprEssence
To assess how network interactions may be predictive of health deterioration we utilised the ExprEssence method of network condensation3 to reduce a STRING protein-protein interaction network to smaller, disease-relevant subnetworks. ExprEssence works by mapping expression values to a previously constructed gene network such as those provided by STRING.27 Each gene in the network is annotated with two expression values (eg, one from averaging all control expression values and one from averaging all case expression values) and then a ‘link score’ is calculated for each of that gene’s edges in the network by aggregating the differences between the control and case values of both genes at either side of the edge. The network can then be condensed by setting a threshold for this link score which an edge needs to meet in order to remain in the network or by taking the top-N edges with the highest link scores.3 We used the STRING database, selecting only physical interactions and those with a quality score ≥999 so as to create an interaction network of only those interactions with the highest confidence (a total of 9691 interactions). After calculating the link scores of the samples (making the relevant case vs control comparison for each dataset), we condensed the network down to the default of 50 edges with the highest link scores, calculated the average link score for each subnetwork and used this average to rank them. These 50 edges produced subnetworks and (up to) the first 5 were taken as the ExprEssence features.
Random survival forest
Following the SASKit-ML pipeline as closely as possible, where features could not be calculated, they were left out of the analysis (overview in online supplemental table S2). Integrated features were used to train a random survival forest model, using the randomForestSRC package of R, with evaluation carried out using the survcomp package to calculate the Integrated Brier Score (IBS) and C-index using 10-fold cross-validation, following the example in Dereli et al.28 An IBS of 0 indicates perfect prediction, and 0.25 indicates the model is no better than guessing, while a C-index of 1 indicates perfect prediction and 0.5 indicates the model is no better than guessing. During the cross-validation, features were integrated independently for each fold with the reported features for each model being those calculable by using the complete dataset. For the HMP and RA-MAP datasets, the variable importance (VIMP) score of each feature is also reported.
Supplemental material
Patient and public involvement
None.
Results
HMP and RA-MAP
Based on the 10-fold cross-validation, the T2D predictive model had a modest predictive ability with a mean C-index of 0.632 and a mean IBS of 0.107. The RA model was similar, with a mean C-index of 0.654 and a mean IBS of 0.150. These values indicate a modest ability to predict health deterioration in patients with potential pre-diabetes and RA. Finally, the IBD model was no better than guessing when assessed with C-index (mean of 0.507) or by IBS (mean of 0.238). It should be noted that variation across the 10-folds was high for both the RA and T2D models, likely due to the small sample sizes of these datasets (figure 1).
Out of the 25 features used in the T2D model, 6 had positive, non-negligible VIMP scores. Age proved to be the most predictive, alongside the transcriptomic GSVA activation score for ‘Arrhythmogenic Right Ventricular Cardiomyopathy ARVS’ and hsCRP (figure 2A) and the proteomic ExprEssence subnetwork 2 (HRG–PLG) (figure 3A).
Out of the 13 features used in the IBD model, 4 had positive, non-negligible VIMP scores, the transcriptomic ExprEssence subnetworks 4 (CSF3R–CSF3) and 3 (CXCL1, CXCL2, CXCL5, CXCR2) (figure 3B), and the transcriptomic GSVA activation score for the KEGG pathway ‘Phenylalanine Metabolism’ being the most important (figure 2B).
Out of the 20 features used in the RA model, 4 had positive VIMP scores. The most predictive by far was the transcriptomic GSVA activation score for the KEGG pathway ‘Primary Immunodeficiency’, followed by neutrophil-lymphocyte ratio and age (figure 2C).
The Cancer Genome Atlas
Across TCGA data, the only consistent data available appropriate for the SASKit-ML pipeline were the endpoint and transcriptomics data, and as such, to allow comparability between the different cancers, only transcriptomic features, age and sex were used.
Most TCGA produced models with an average C-index between 0.55 and 0.7 (figure 4A) and an average IBS below 0.2 (figure 4B), however, some models, particularly those for cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), pancreatic adenocarcinoma (PAAD) and rectum adenocarcinoma (READ) saw substantial variation in the 10-fold cross-validation of one or both statistics, suggesting possible overfitting in these datasets in some folds.
Discussion
Presented in this paper are the results from applying the SASKit-ML biomarker learning pipeline to several public datasets: the T2D and IBD branches of the HMP, the RA-MAP and TCGA. The 19 total datasets encompass pre-diabetes, IBD, RA and cancer. The principal challenge in adapting the SASKit-ML pipeline to these datasets was in drawing parity between the available features, and those available in the SASKit cohort that motivated the prepublished data analysis plan, with the HMP-T2D dataset matching most closely, followed by RA-MAP and the HMP-IBD. TCGA data were the least comparable to SASKit. Attempting to adhere as closely as possible to the published SASKit-ML pipeline is both the main strength and limitation of this study. We only deviated from the published protocol when it was demanded by the differences in study cohorts to the SASKit cohort itself, mostly due to differing data types and sample sizes. This means that, despite not having separate testing data to validate our models and having to rely on 10-fold cross-validation of the training data, we have not risked the cherry-picking of features as a source of overfitting. However, close adherence to the published protocol also means that we were unable to be flexible in adapting the pipeline to the new datasets, for instance, by adding in known clinically relevant features for each disease, which could have possibly improved the models.
The results show that while the SASKit-ML can produce informative features, they are likely to require further selection (eg, by minimal depth thresholding29). Further, as acknowledged in the study protocol, such a strict data integration approach carries the risk of significant information loss, potentially sacrificing performance in order to improve human interpretability. With that said, none of the analysed datasets were perfectly comparable to the SASKit cohort, lacking some or all of the clinical markers or lacking suitable controls (ie, healthy controls as opposed to healthy tissue adjacent to a tumour or patients presenting with disease elements insufficient to receive a formal diagnosis). As such it can be expected that the pipeline will perform at least as well but likely better on the SASKit cohort, suggesting that our expected model accuracy of around 75% based on the power analysis in the study protocol1 is within reach. Of concern is that, for all models, the cross-validation results tended to be highly variable indicating a propensity towards overfitting, however, this is perhaps unavoidable given the low sample sizes and complex endpoints of the studies.
Focusing on the three datasets most similar to the SASKit cohort, HMP-T2D, HMP-IBD and RA-MAP, a few observations can be made. First, while none of the models produced were outstanding, both the HMP-T2D and RA-MAP models were reasonably predictive. While these two models did not reach the predictive power aimed for in the SASKit project, this may be accounted for by the differences in the study designs and, at least in the case of HMP-T2D, the less clearly defined controls. Importantly, the HMP-IBD endpoint was based on highly variable, subjective assessments of disease progression and quality of life, whereas the HMP-T2D endpoint was based purely on measured laboratory values while the RA-MAP endpoint was based on a mix of symptom-related and quality-of-life-related data types. As such, it is perhaps not surprising that the HMP-T2D model was both more accurate than the other two, but also that the cross-validation results were less variable (figure 1).
The features identified for all three of these datasets broadly made biological sense. Age and hsCRP were both highly predictive of onset of pre-diabetes in the HMP-T2D study, both of which are known predictors of pre-diabetes and known to correlate with each other.30 Proteomic ExprEssence Subnetwork 2 (PLG–HRG) also ties into current understanding of T2D, representing modulation of the plasminogen/histidine rich glycoprotein relationship, implicated in aberrant blood clotting (for which T2D is a risk factor).31
The top feature of the HMP-IBD model was the transcriptomic ExprEssence Subnetwork 4 linkscore between CSF3 and CSF3R, with increased expression shown in the patients with IBD for both genes, which is a commonly observed expression change in IBD studies.32 33 Further, these genes have both been identified as biomarkers of response to treatment in ulcerative colitis.34 The transcriptomic ExprEssence Subnetwork 3 genes (CXCL1, CXCL2, CXCL5 and CXCR2) are markers of general inflammation but also appear to be specifically modulated in ulcerative colitis patients.35 All six genes of these two subnetworks, and thus their interactions, are upregulated in diseased patients compared with the healthy controls, potentially indicating an upregulation of specific immune genes in both ulcerative colitis and Crohn’s disease. The one other positively predictive variable for IBD was the GSVA activation score for the KEGG pathway ‘Phenylalanine Metabolism’, a process known to be altered in IBD in both mice and humans.36 37
Of interest to SASKit as a project interested in the inflammation-related comorbidities of many diseases, including cancer and stroke, it has been shown recently that early elevation of poststroke serum CSF3 correlates with stroke severity and worse 3-month functional outcome38 while the therapeutic application of CSF3 may alleviate cerebral ischaemic injury.39 Further, It has also been shown that the CXCL8-CXCR2 chemotactic axis is upregulated in the peripheral blood of ischaemic stroke patients, with higher serum CXCL8 levels in patients with unfavourable functional outcomes after 3 months.40
Finally, the top features for the RA-MAP study included activation of the transcriptomic GSVA score for the KEGG pathway ‘Primary Immunodeficiency’, a set of conditions closely connected to RA,41 with expression of genes from this pathway having been shown to be good discriminators of RA from osteoarthritis.42 Additionally, neutrophil-lymphocyte ratio and age were both top features for predicting the worsening of RA, both of which were established biomarkers of RA prognosis.43 44 Further, neutrophil-lymphocyte ratio has been shown by meta-analysis of ~2500 patients and controls to be a strong indicator of the presence of RA.45
While most of the top features had clear connections to the disease context they were selected for, this was not always the case. The GSVA score for the arrhythmogenic right ventricular cardiomyopathy KEGG pathway for instance was the second most important feature for the HMP-T2D model. This KEGG pathway contains the genes SGCA, SGCB, SGCD and SGCG which code proteins of the sarcoglycan complex. SGCG has been shown to be associated with T2D in a Punjabi Sikh population in India46 while SGCZ (related but not present in this KEGG term) has been associated with pre-diabetes status change.47
Given the overall poor to middling performance of the models, particularly those in the datasets closest to the SASKit cohort, the question is raised as to whether the SASKit-ML pipeline could be improved or adapted not just to these datasets, but to the SASKit cohort itself. One notable oversight in the SASKit-ML pipeline is that it does not include a feature elimination step to remove uninformative features (which were included in the models presented here, but not visualized in the VIMP plots). One possibility to improve the resultant SASKit models would thus be to apply a feature elimination step such as minimal depth thresholding.48 This could potentially improve results where the pipeline generated many unimportant features, for instance, in the HMP-IBD and RA-MAP cohorts where only four of the features had positive VIMP scores (figure 2). Further, this study raised the potential issue of low signal between patients and controls. It could be that a high degree of difference between patient and control omics is required for the features to be predictive, and in the case of the HMP cohorts here, the controls were not neatly defined and signal was low. If this problem occurs in the SASKit cohort as well, it could be worth considering alternate comparisons for the feature generation in the exploratory analyses, for instance, by comparing the omics of patients who we retrospectively know to suffer health deterioration to those who do not.
Overall, this analysis has established that the published SASKit-ML pipeline does have some generalisability to other diseases, but that the resultant models can be highly variable and rarely highly predictive. With that said, the features produced, especially those that demonstrated some predictive power, were biologically plausible, broadly explainable and able to provide some insight into potential mechanisms of the studied diseases.
Data availability statement
Data are available in a public, open access repository. Data are available on reasonable request. The data used in this study come from the publicly available datasets of the HMP, RA-MAP and TCGA. R scripts and preprocessed data to replicate the study are provided at https://github.com/danrpalmer/saskit-ml.
Ethics statements
Patient consent for publication
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
Contributors DP and GF conceived of the study. DP performed the analyses and wrote the manuscript. DP, LH, HME, UW, AK and GF contributed to the analysis plan, discussion and review of the manuscript. The guarantor of the study is GF; who accepts full responsibility for the finished work and the conduct of the study, had access to the data and controlled the decision to publish. ChatGPT was used to help streamline code and to assist with shortening the manuscript.
Funding We acknowledge the financial support by the Federal Ministry of Education and Research (BMBF) of Germany for the SASKit study (FKZ 01ZX1903A).
Competing interests UW reports grants and personal fees from Merz Pharma, personal fees from Amarin, personal fees from Bristol-Myers Squibb, personal fees from Canon Medical Systems, personal fees from Daiichi Sankyo, personal fees from Ipsen Pharma, personal fees from Pfizer, personal fees from Thieme and personal fees from Elsevier Press, all outside the submitted work.
Patient and public involvement Patients and/or the public were not involved in the design, or conduct, or reporting, or dissemination plans of this research.
Provenance and peer review Not commissioned; externally peer reviewed.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.