Correspondence to Dr Georg Fuellen; [email protected]
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.
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.
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).
Enlarge this image.Figure 1. Prediction accuracy of health deterioration using the C-index (A) and Integrated Brier Score (B) for the T2D and IBD datasets of the HMP and the RA dataset of RA-MAP. The points correspond to the results of different folds in the 10-fold cross-validation. HMP, Human Microbiome Project; IBD, inflammatory bowel disease; IBS, Integrated Brier Score; RA-MAP, Rheumatoid Arthritis Map; T2D, type 2 diabetes.
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).
Enlarge this image.Figure 2. Positive VIMP scores for the features of the SASKit model constructed for the T2D (A) and IBD (B) datasets of the HMP and the RA (C) dataset of RA-MAP. HMP, Human Microbiome Project; hs-CRP, high-sensitive C reactive protein; IBD, inflammatory bowel disease; RA, rheumatoid arthritis; T2D, type 2 diabetes; VIMP, variable importance.
Enlarge this image.Figure 3. ExprEssence subnetworks with positive VIMP scores for the T2D (A) and IBD (B) datasets of the HMP. On the left side of each box, expression in healthy samples is visualised and on the right side of each box, expression in disease samples is visualised; the boxes are connected by the corresponding link-scores. In both figures, red boxes indicate high expression relative to average gene expression while green indicates low expression relative to average gene expression, with white indicating an average expression level. Note that in (A) both sides of the boxes are saturated red as both genes were highly expressed, the (weak) linkscore is positive as expression of both HRG and PLG was slightly elevated in the patients with pre-diabetes relative to the controls. HMP, Human Microbiome Project; IBD, inflammatory bowel disease; T2D, type 2 diabetes; VIMP, variable importance.
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.
Enlarge this image.Figure 4. Prediction accuracy of health deterioration using the C-index (A) and Integrated Brier Score (IBS) (B) for the TCGA datasets over 10-fold cross-validation. TCGA, The Cancer Genome Atlas.
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
Not applicable.
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.
1 Henze L, Walter U, Murua Escobar H, et al. Towards biomarkers for outcomes after pancreatic ductal adenocarcinoma and ischaemic stroke, with focus on (co)-morbidity and ageing/cellular senescence (SASKit): protocol for a prospective cohort study. BMJ Open 2020; 10: e039560. doi:10.1136/bmjopen-2020-039560
2 Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013; 14: 7. doi:10.1186/1471-2105-14-7
3 Warsow G, Greber B, Falk SSI, et al. ExprEssence--revealing the essence of differential experimental data in the context of an interaction/regulation net-work. BMC Syst Biol 2010; 4: 164. doi:10.1186/1752-0509-4-164
4 Lin DY. Cox regression analysis of multivariate failure time data: the marginal approach. Stat Med 1994; 13: 2233–47. doi:10.1002/sim.4780132105
5 Ishwaran H, Kogalur UB, Blackstone EH, et al. Random survival forests. Ann Appl Stat 2008; 2: 841–60. doi:10.1214/08-AOAS169
6 López de Maturana E, Alonso L, Alarcón P, et al. Challenges in the Integration of Omics and Non-Omics Data. Genes (Basel) 2019; 10: 238. doi:10.3390/genes10030238
7 Geistlinger L, Csaba G, Santarelli M, et al. Toward a gold standard for benchmarking gene set enrichment analysis. Brief Bioinform 2021; 22: 545–56. doi:10.1093/bib/bbz158
8 Herrmann M, Probst P, Hornung R, et al. Large-scale benchmark study of survival prediction methods using multi-omics data. Brief Bioinform 2021; 22: bbaa167. doi:10.1093/bib/bbaa167
9 Si Z, Sun L, Wang X. Evidence and perspectives of cell senescence in neurodegenerative diseases. Biomed Pharmacother 2021; 137: 111327. doi:10.1016/j.biopha.2021.111327
10 Mehdizadeh M, Aguilar M, Thorin E, et al. The role of cellular senescence in cardiac disease: basic biology and clinical relevance. Nat Rev Cardiol 2022; 19: 250–64. doi:10.1038/s41569-021-00624-2
11 Schafer MJ, Zhang X, Kumar A, et al. The senescence-associated secretome as an indicator of age and medical risk. JCI Insight 2020; 5: e133668. doi:10.1172/jci.insight.133668
12 Martyanov V, Whitfield ML, Varga J. Senescence Signature in Skin Biopsies From Systemic Sclerosis Patients Treated With Senolytic Therapy: Potential Predictor of Clinical Response? Arthritis Rheumatol 2019; 71: 1766–7. doi:10.1002/art.40934
13 Creasy HH, Felix V, Aluvathingal J, et al. HMPDACC: a Human Microbiome Project Multi-omic data resource. Nucleic Acids Res 2021; 49: D734–42. doi:10.1093/nar/gkaa996
14 Proctor LM. The Integrative Human Microbiome Project. Nature New Biol 2019; 569: 641–8. doi:10.1038/s41586-019-1238-8
15 Narasimhan A, Flores RR, Robbins PD, et al. Role of Cellular Senescence in Type II Diabetes. Endocrinology 2021; 162: bqab136. doi:10.1210/endocr/bqab136
16 Choi EL, Taheri N, Chandra A, et al. Cellular Senescence, Inflammation, and Cancer in the Gastrointestinal Tract. Int J Mol Sci 2023; 24: 9810. doi:10.3390/ijms24129810
17 Singh V, Kalliolias GD, Ostaszewski M, et al. RA-map: building a state-of-the-art interactive knowledge base for rheumatoid arthritis. Database (Oxford) 2020; 2020: baaa017. doi:10.1093/database/baaa017
18 Isaacs JD, Brockbank S, Pedersen AW, et al. RA-MAP, molecular immunological landscapes in early rheumatoid arthritis and healthy vaccine recipients. Sci Data 2022; 9: 196. doi:10.1038/s41597-022-01264-y
19 TCGA Research Network. Data from: The cancer genome atlas,TCGA. 2019. Available: https://www.cancer.gov/ccg/research/genome-sequencing/tcga
20 American Diabetes Association. Diagnosis and classification of diabetes mellitus. Diabetes Care 2013; 36 Suppl 1: S67–74. doi:10.2337/dc13-S067
21 Allegretti JR, Barnes EL, Stevens B, et al. Predictors of Clinical Response and Remission at 1 Year Among a Multicenter Cohort of Patients with Inflammatory Bowel Disease Treated with Vedolizumab. Dig Dis Sci 2017; 62: 1590–6. doi:10.1007/s10620-017-4549-3
22 Barnes EL, Liew C-C, Chao S, et al. Use of blood based biomarkers in the evaluation of Crohn’s disease and ulcerative colitis. World J Gastrointest Endosc 2015; 7: 1233–7. doi:10.4253/wjge.v7.i17.1233
23 Edwards CJ, Kiely P, Arthanari S, et al. Predicting disease progression and poor outcomes in patients with moderately active rheumatoid arthritis: a systematic review. Rheumatol Adv Pract 2019; 3: rkz002. doi:10.1093/rap/rkz002
24 Liu J, Lichtenberg T, Hoadley KA, et al. An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics. Cell 2018; 173: 400–16. doi:10.1016/j.cell.2018.02.052
25 Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005; 102: 15545–50. doi:10.1073/pnas.0506580102
26 Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000; 28: 27–30. doi:10.1093/nar/28.1.27
27 Szklarczyk D, Kirsch R, Koutrouli M, et al. The STRING database in 2023: protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res 2023; 51: D638–46. doi:10.1093/nar/gkac1000
28 Dereli O, Oğuz C, Gönen M. Path2Surv: Pathway/gene set-based survival analysis using multiple kernel learning. Bioinformatics 2019; 35: 5137–45. doi:10.1093/bioinformatics/btz446
29 Ishwaran H, Kogalur UB, Gorodeski EZ, et al. High-dimensional variable selection for survival data. J Am Stat Assoc 2010; 105: 205–17. doi:10.1198/jasa.2009.tm08622
30 Ghule A, Kamble TK, Talwar D, et al. Association of Serum High Sensitivity C-Reactive Protein With Pre-diabetes in Rural Population: A Two-Year Cross-Sectional Study. Cureus 2021; 13: e19088. doi:10.7759/cureus.19088
31 Szułdrzyński K, Jankowski M, Potaczek DP, et al. Plasma Fibrin Clot Properties as Determinants of Bleeding Time in Human Subjects: Association with Histidine-Rich Glycoprotein. Dis Markers 2020; 2020: 7190828. doi:10.1155/2020/7190828
32 Hotte NSC, Salim SY, Tso RH, et al. Patients with inflammatory bowel disease exhibit dysregulated responses to microbial DNA. PLoS One 2012; 7: e37932. doi:10.1371/journal.pone.0037932
33 Liu Y, Duan Y, Li Y. Integrated Gene Expression Profiling Analysis Reveals Probable Molecular Mechanism and Candidate Biomarker in Anti-TNFα Non-Response IBD Patients. J Inflamm Res 2020; 13: 81–95. doi:10.2147/JIR.S236262
34 Haberman Y. Tissue-based Gene Expression as Potential Biomarkers for IBD Course. Inflamm Bowel Dis 2020; 26: 1485–9. doi:10.1093/ibd/izaa217
35 Shi L, Han X, Li J-X, et al. Identification of differentially expressed genes in ulcerative colitis and verification in a colitis mouse model by bioinformatics analyses. World J Gastroenterol 2020; 26: 5983–96. doi:10.3748/wjg.v26.i39.5983
36 Wikoff WR, Anfora AT, Liu J, et al. Metabolomics analysis reveals large effects of gut microflora on mammalian blood metabolites. Proc Natl Acad Sci U S A 2009; 106: 3698–703. doi:10.1073/pnas.0812874106
37 Bosch S, Struys EA, van Gaal N, et al. Fecal Amino Acid Analysis Can Discriminate De Novo Treatment-Naïve Pediatric Inflammatory Bowel Disease From Controls. J Pediatr Gastroenterol Nutr 2018; 66: 773–8. doi:10.1097/MPG.0000000000001812
38 Wall A, Anger O, Jood K, et al. Circulating granulocyte colony-stimulating factor and functional outcome after ischemic stroke: an observational study. Neurol Res 2021; 43: 1013–22. doi:10.1080/01616412.2021.1948766
39 Patel AMR, Apaijai N, Chattipakorn N, et al. The Protective and Reparative Role of Colony-Stimulating Factors in the Brain with Cerebral Ischemia/Reperfusion Injury. Neuroendocrinology 2021; 111: 1029–65. doi:10.1159/000512367
40 He Q, Shi X, Zhou B, et al. Interleukin 8 (CXCL8)-CXC chemokine receptor 2 (CXCR2) axis contributes to MiR-4437-associated recruitment of granulocytes and natural killer cells in ischemic stroke. Mol Immunol 2018; 101: 440–9. doi:10.1016/j.molimm.2018.08.002
41 Sogkas G, Witte T. The link between rheumatic disorders and inborn errors of immunity. EBioMedicine 2023; 90: 104501. doi:10.1016/j.ebiom.2023.104501
42 Zhang R, Yang X, Wang J, et al. Identification of potential biomarkers for differential diagnosis between rheumatoid arthritis and osteoarthritis via integrative genome‑wide gene expression profiling analysis. Mol Med Rep 2019; 19: 30–40. doi:10.3892/mmr.2018.9677
43 Mercan R, Bitik B, Tufan A, et al. The Association Between Neutrophil/Lymphocyte Ratio and Disease Activity in Rheumatoid Arthritis and Ankylosing Spondylitis. J Clin Lab Anal 2016; 30: 597–601. doi:10.1002/jcla.21908
44 Innala L, Berglin E, Möller B, et al. Age at onset determines severity and choice of treatment in early rheumatoid arthritis: a prospective study. Arthritis Res Ther 2014; 16: R94. doi:10.1186/ar4540
45 Erre GL, Paliogiannis P, Castagna F, et al. Meta-analysis of neutrophil-to-lymphocyte and platelet-to-lymphocyte ratio in rheumatoid arthritis. Eur J Clin Invest 2019; 49: e13037. doi:10.1111/eci.13037
46 Saxena R, Saleheen D, Been LF, et al. Genome-wide association study identifies a novel locus contributing to type 2 diabetes susceptibility in Sikhs of Punjabi origin from India. Diabetes 2013; 62: 1746–55. doi:10.2337/db12-1077
47 Liu T, Li H, Conley YP, et al. A Genome-Wide Association Study of Prediabetes Status Change. Front Endocrinol (Lausanne) 2022; 13: 881633. doi:10.3389/fendo.2022.881633
48 Gilhodes J, Zemmour C, Ajana S, et al. Comparison of variable selection methods for high-dimensional survival data with competing events. Comput Biol Med 2017; 91: 159–67. doi:10.1016/j.compbiomed.2017.10.021
© 2024 Author(s) (or their employer(s)) 2024. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ. http://creativecommons.org/licenses/by-nc/4.0/ 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/ . Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Details
1 Institute for Biostatistics and Informatics in Medicine and Ageing Research, Rostock University Medical Center, Rostock, Germany
2 Department of Medicine, Clinic III, Hematology, Oncology, Palliative Medicine, Rostock University Medical Center, Rostock, Germany; Department of Internal Medicine II - Hematology, Oncology and Palliative Medicine, Asklepios Hospital Group Harz Mountains, Goslar, Germany
3 Department of Medicine, Clinic III, Hematology, Oncology, Palliative Medicine, Rostock University Medical Center, Rostock, Germany
4 Department of Neurology, Rostock University Medical Center, Rostock, Germany