Using serum metabolomics analysis to predict sub-clinical atherosclerosis in patients with SLE

Patient cohort

Serum samples were collected from 80 non-fasting patients with SLE attending a rheumatology clinic at University College London Hospital (UCLH) and fulfilling the American College of Rheumatology (ACR) classification criteria for lupus (1997)[18]. Patients had no previous history of CVD (defined as coronary artery disease, stroke, or myocardial infarction with confirmatory evidence from blood tests and/or imaging) and were scanned by vascular ultrasound between 2011-2013 [19].

Carotid and femoral ultrasound scans were performed to identify the presence of plaques. Demographic and clinical information were recorded at the time of scan/blood sampling, including sex, age, ethnicity, mean blood pressure (BP, average between systolic and diastolic BP), treatment (including hydroxychloroquine (HCQ), statins, ACE inhibitor, immunosuppressives, rituximab, prednisolone (and dose), and aspirin), and disease activity assessed by the global British Isles Lupus Assessment Group (BILAG) index [20] (Table 1). In total, four patients were not on any treatment at the time of the scan. All patients gave informed written consent and the study was approved by the combined UCL/UCLH Research Ethics Committee (Reference 06/Q0505/79).

Plaque detection

Scans were performed by an experienced vascular scientist using strict scanning protocols and established equipment settings (Griffin et al., 2009; Nicolaides et al., 2010). Vascular ultrasound scans of the common carotid artery, carotid bulb, carotid bifurcation, common femoral artery and femoral bifurcation were performed bilaterally using the Philips IU22 ultrasound computer and the L9-3 MHz probe. IMT measurements were performed using QLAB Advanced Quantification Software® version 7.1 (Philips Ultrasound, Bothell, USA). Plaque was defined as “a focal thickening >1.2 mm that encroaches into the arterial lumen as measured from the media-adventitia interface to the lumen interface” [21]. Patients having at least one region fulfilling this description were included in the group with plaque (SLE-P).

Ultrasound images were stored as DICOM files and analysed using Carotid Plaque Texture Analysis software for ultrasonic arterial wall and atherosclerotic plaque measurements (LifeQ Medical Ltd – www.lifeqmedical.com). Total plaque area (TPA) was defined as the sum of the cross-sectional areas of all plaques seen in longitudinal images (plaque area in mm2). Grey Scale Median (GSM): Plaque echogenicity was expressed numerically by GSM value [22], a measure of plaque stability and lipid content. Lower GSM values signify more echolucent plaque associated with lipid and inflammatory cell content; whereas plaque with high GSM scores are more stable and less echolucent [23]. Images were normalized using linear scaling with two reference points blood (grey scale=0) and adventitia (grey scale=190).

Serum metabolomics analysis

Measures of 228 serum biomarkers were acquired with an established nuclear magnetic resonance (NMR)-spectroscopy platform (Nightingale Health)[24, 25]. These included both absolute concentrations (mmol/L), ratios, and percentages (%) of lipoprotein composition. Serum lipids measured included apolipoproteins (Apo) and VLDL, LDL, intermediate density lipoprotein (IDL) and HDL particles of different sizes ranging from chylomicrons and extremely large (XXL), very large (XL), large (L), medium (M), small (S) and very small (XS). Lipids within each lipoprotein subclass included –total lipid (L), phospholipids (PL), total cholesterol (C), cholesterol esters (CE), free cholesterol (FC) and triglycerides (TG). Distribution of these lipids was expressed as a ratio or percentage of total lipid content for each lipoprotein subclass (for list of metabolites see Supplementary Table 1).

Logistic regression

To assess the association of plaque with NMR metabolomic biomarker data, univariate logistic regressions were performed for each individual serum metabolite, adjusted for ethnicity, age, mean blood pressure, global BILAG score, and treatments at the time of the scan. Results were visualized in a forest plot using R package forestplotNMR [26].

Predictive models

The data analysis pipeline used in shown in Figure 1. RStudio (The R Foundation, Vienna, Austria) (R Core Team, 2019) and Orange 3.24.1 (Bioinformatics Lab, University of Ljubljana, Slovenia) [27] were used for machine learning analysis as we have shown previously [28]. Five different supervised learning algorithms were implemented: support vector machine (SVM), logistic regression with and without interactions, decision trees, and random forest classification. The outcome of the learning algorithms was to predict whether an SLE patient is likely to develop plaque. Predictive models were generated from metabolite concentrations at the time of the scan.

• Download figure
• Open in new tab

Figure 1. Data analysis workflow

Flow chart depicting the data cleaning and processing steps taken prior data analysis using machine learning algorithms. Abbreviations: SVM – support vector machine | P – plaque | NP – no plaque, BILAG – British Isles Lupus Assessment Group disease activity score.

Missing data: Features with >10% missing data were excluded (pyruvate). Remaining missing values (n=5) were imputed using k-nearest neighbors with k = 5.

Homology reduction: To reduce homology, if two features had a correlation co-efficient >0.95 then the feature with the greatest mean absolute correlation with the remaining features was removed (Supplementary Data File 1). This left 124 metabolites.

Data scaling: Metabolite concentrations were centered on the mean and scaled to the standard deviation.

Predictors: The independent variables included in the models were the homology reduced dataset (124 metabolites), as well as the cohort information (age, ethnicity, mean blood pressure, disease duration, BILAG, treatments at the time of first scan). Sex was not considered as all participants were female. Full lists of the predictors contributing to each model are included in Supplementary Data File 2.

Support vector machine (SVM): A supervised classification method which creates a hyperplane to optimally separate data into two classes [29]. As this data set was not linearly separable, the radial basis function kernel was used. Values for C, epsilon, and gamma were tuned using the R Package e1071 [30]. The parameters were set to C = 4.0, epsilon = 0.1, gamma = 0.01.

Decision tree: A supervised machine learning method which classifies incidents according to their features. Decision trees are built using forms of impurity measures, such as information gain and entropy [31]. To prevent overfitting, decision trees were limited to a depth of 4 and subsets of 5 or less were not split further.

Random Forest (RF): A statistical classifier (machine-learning algorithm) that assigns observations into classes by creating a set of decision trees, or ‘forest’. Importance of each feature was quantified by the Gini index, which represents the total variance across the two classes, the purity of each node and the quality of each split. The values for mtry and ntree were tuned using the R Package [32]. The parameters were set to mtry=10, ntree=10,000.

Logistic regression with/without interactions (LR/LR+I): The least absolute shrinkage and selection operator (lasso) method uses the absolute value of the co-efficient as a penalty to shrink less important features to zero. The strength of shrinkage is determined by tuning the regularization variable lambda (λ). Ln(λ) was optimized using the R package [32] and set to lambda=0.059 and lambda=0.065 for logistic regression with and without interaction respectively.

Model performance: Ten-fold cross-validation was used to evaluate model performance. Validation was performed in Orange for the decision tree and in R for all other models, using a balanced splitting. The following performance metrics were calculated from the confusion matrices: (i) F1 score - a weighted average of precision (positive predictive value) and recall (sensitivity), (ii) specificity - the true negative rate, and (iii) classification accuracy (CA) - the proportion of correctly classified cases.

Partial Least Squares Discriminant Analysis (sPLS-DA)

A sPLS-DA is a supervised clustering machine-learning approach that combines parameter selection and classification into one operation. This analysis was performed using the features selected for the lasso LR+I model (Table 3).

Statistical testing

Statistical tests were performed in Microsoft Excel and GraphPad Prism version 8.3.0 for Windows (GraphPad Software, San Diego, USA). Data was assessed for normality and analysed with parametric or non-parametric tests as appropriate. Details of statistical tests are given in the figure legends. P-values < 0.05 were considered statistically significant.

Serum metabolites can predict the presence of pre-clinical atherosclerotic plaque in SLE patients

Metabolites were quantified in serum from SLE patients classed as having plaque(s) (SLE-P) or not having plaque(s) (SLE-NP) using vascular ultrasound (Table 1 and Supplementary Table 1). Several analysis strategies were applied to identify metabolites associated with sub-clinical plaque development (Figure 1). All models were adjusted for ethnicity, age, mean blood pressure (BP), disease duration, disease activity (global BILAG score), and treatment at the time of the scan. Firstly, univariate logistic regressions identified four metabolites which differentiated between SLE-P and SLE-NP patients; Leucine, M-VLDL-FC_%, L-VLDL-FC_% and XXL-VLDL-PL_% which were all increased in serum from SLE-P compared to SLE-NP patients (Figure 2A, Supplementary Data File 3 and 4). Of note, these metabolites can be significantly affected by treatment with statins [33, 34] and pro-protein convertase subtilisin/kexin type 9 (PCSK9) inhibitors [35], or body mass index (BMI) [36], suggesting that abnormal serum lipid metabolite profiles in SLE-P patients, could be modified using available therapies or interventions (Figure 2A asterisks).

• Download figure
• Open in new tab

Figure 2. Identification of important metabolites separating patients with SLE-P from SLE-NP

(A) Forest plot depicting statistically significant individual logistic regression results of metabolites in SLE-P vs SLE-NP. Results given in Odds Ratio (95% confidence intervals). Logistic regressions were adjusted for age, ethnicity, mean blood pressure, global BILAG, disease duration, and treatments at the time of scan. Coloured asterisks denote metabolite has previously been shown to be modified by statins, PCSK9 inhibitors or BMI [33-36]. (B) Best performing models were determined based on performance statistics (Table 2). The top ten features of these LR and RF models are listed to identify common features. (C-D) Metabolites which were significant in the individual logistic regression and featured in the top ten of the LR and RF models were further analysed using (C) bar charts showing mean, and p value (see Supplementary Figure 2A) and (D) ROC plots (see Supplementary Figure 2B). (E) Metabolites which featured in the top ten of the LR and RF model were combined and analysed in a single ROC plot (L-VLDL-FC%, Glycine, XXL-VLDL-PL%, Tyrosine). (F) Beta coefficients from the logistic regressions with interactions are plotted for the SLE-P vs SLE-NP analysis. The sign indicates the direction of the effect of the predictor; a positive sign indicates an increased likelihood of a SLE-P prediction, while a negative sign indicates an increased likelihood of a SLE-NP prediction.

Next, five supervised ML models were developed and validated to predict the presence of plaque in SLE patients; LR, LR+I, SVM, RF and Decision tree. Since many of the metabolites measured were biologically interdependent, and therefore highly correlated, homology reduction was applied (see methods and Supplementary Data File 1). The models were built using the homology reduced dataset (124 metabolites) and patient information (age, ethnicity, mean BP, disease duration, global BILAG index, treatments at the time of first scan). Sex was not considered as all participants were female (Supplementary Data File 2 for full lists of the predictors contributing to each model).

The top three models, according to classification accuracy, specificity, and F1 scores, were LR, LR+I, and RF (Table 2). Performance metrics were based on predictions of the models, summarised in confusion matrices (Supplementary Figure 1). LR and LR+I had a similar performance, correctly classifying 75% and 80% of SLE-P patients respectively. The RF model had the best specificity, identifying 45 out of 50 (90%) SLE-NP patients correctly. The three models were further investigated to identify the top features in predicting plaque formation in SLE (Figure 2B-E). Four metabolites (XXL-VLDL-PL_%, L-VLDL-FC_%, glycine and tyrosine) and patient age were identified in both LR and RF models as important predictors for SLE-P (Figure 2B). Metabolites which were significant in individual logistic regressions and also featured in the top ten metabolites of the LR and RF model, were further investigated for differences between SLE-P and SLE-NP patients (Figure 2C and Supplementary Figure 2A). Receiver operating characteristic curve (ROC) of these metabolites showed an area under the curve (AUC) of 0.6810 (XXL-VLDL-PL %), 0.7523 (L-VLDL-FC %), 0.7337 (glycine), and 0.6300 (tyrosine) (Figure 2D and Supplementary Figure 2B), demonstrating a potential predictive ability of these metabolites independent of other clinical features and metabolites. When the four shared metabolites between LR and RF models were combined, the ROC analysis showed an AUC of 0.7667 (Figure 2E), demonstrating a stronger potential predictive ability of the metabolites grouped together.

The LR+I model included interaction of each metabolite with all other metabolites and clinical features (assessing over 15,000 possible features) and was therefore assessed separately from the other models. This model performed best overall, with a classification accuracy of 0.800 (Table 2). Using the lasso method of shrinkage and selection, only 35 interactions were identified as important and given a beta coefficient to describe the effect size and direction on the model (Figure 2F, Table 3). Features with larger beta coefficients had the greatest effect on the classification. Notably, L-VLDL-FC_%: age and IDL-TG_%:age, Black/Carribean ethnicity:hydroxycloroquine treatment and sphingomyelin:leucine were significantly associated with plaque; while S-LDL-CE:rituximab treatment, S-HDL-FC:XXL-VLDL-CE_%, glycine:histidine and statins:ACE inhibitor treatment were associated with no-plaque.

Metabolite interactions can differentiate between SLE-P and SLE-NP patients

Using the metabolite interaction features that were selected for the LR+I model (Table 3, Figure 2F), a sPLS-DA was performed to rank and validate the metabolite features by their distribution in patients with SLE-P and SLE-NP. By assessing the overall estimation error rate in ten-fold cross-validation, models with 4 components and a subset of 28 metabolite features were chosen for optimal model performance (Figure 3A). This analysis identified a significant separation between SLE-P and SLE-NP patients by plotting principal component (PC)-2 against PC-1 (Figure 3B). The 28 selected metabolite interaction features were ranked by discriminating capability (Figure 3C, 3D). The two highest weighted features were L-VLDL-FC_%:Age and IDL-TG_%:Age, which were also the highest ranked features in the LR+I model (Figure 2F). L-VLDL-FC_% and age were both also featured as one of the top ten features shared by LR and RF (Figure 2B) and were differentially expressed between SLE-P and SLE-NP patients (Figure 2A, Table 1) suggesting an influential role of both features in SLE-P patients.

Differential metabolites correlated with clinical features of SLE-P patients

Finally, the top ten metabolites from LR, RF, and all metabolite interactions included in the LR+I models were correlated with clinical and plaque features (Figure 4, Supplementary Table 2). Significant correlations include GSM (grey scale median, measure of plaque stability and lipid content) correlated positively with S-HDL-TG, XS-VLDL-TG, XL-HDL-TG:Glycerol and IDL-CE:Glycerol and negatively with LDL-D, His:XS-VLDL-CE_% and Glycine:XS-VLDL-CE_%; Plaque number and plaque thickness negatively correlated with His:XS-VLDL-CE_%; TPA (total plaque area) positively correlated with L-VLDL-TG:Lactate; and disease activity (BILAG) positively correlated with M-VLDL-FC_% and M-VLDL-FC_%:L-LDL-FC_% and negatively with DHA-FA:age. The strongest correlation was between disease duration and Tyrosine:disease duration.

• Download figure
• Open in new tab

Figure 3. Partial least squares discriminate analysis validated metabolites identified by logistic regression with interactions to predict SLE-P

(A) Model optimisation – model with different components and features kept in the analysis were analysed, with each colour representing a different number of components (Comp), number of features kept in the analysis on the x-axis, and the overall error on the y-axis. (B) sPLS-DA plot to validate top hits from the logistic regression with interactions. sPLS-DA is a supervised clustering method which separates SLE-P from SLE-NP. (C) Features included in the sPLS-DA plotted with their factor loading value. (D) Visualisation of the weighting and correlation of each metabolite in component 1 and 2 on the sPLS-DA model.

• Download figure
• Open in new tab

Figure 4. Correlations of significant metabolites with SLE clinical markers

Correlations between metabolites and patient clinical characteristics. Pearson’s product moment correlation coefficients are represented as connecting lines between the clinical characteristic section (grey) and metabolite section (rainbow). Only correlations with p value below 0.05 are shown. Red line=positive correlation and blue line=negative correlation. Width of lines represents the value of correlation coefficients (measured with scale). GSM – grey scale median, BILAG – British Isles Lupus Group Index disease activity score, Plq No – plaque number, TPT – total plaque thickness, TPA – Total Plaque area. (see Supplementary Table 1 for metabolite abbreviations.