A classification modeling approach for determining metabolite signatures in osteoarthritis

Multiple factors can help predict knee osteoarthritis (OA) patients from healthy individuals, including age, sex, and BMI, and possibly metabolite levels. Using plasma from individuals with primary OA undergoing total knee replacement and healthy volunteers, we measured lysophosphatidylcholine (lysoPC) and phosphatidylcholine (PC) analogues by metabolomics. Populations were stratified on demographic factors and lysoPC and PC analogue signatures were determined by univariate receiver-operator curve (AUC) analysis. Using signatures, multivariate classification modeling was performed using various algorithms to select the most consistent method as measured by AUC differences between resampled training and test sets. Lists of metabolites indicative of OA [AUC > 0.5] were identified for each stratum. The signature from males age > 50 years old encompassed the majority of identified metabolites, suggesting lysoPCs and PCs are dominant indicators of OA in older males. Principal component regression with logistic regression was the most consistent multivariate classification algorithm tested. Using this algorithm, classification of older males had fair power to classify OA patients from healthy individuals. Thus, individual levels of lysoPC and PC analogues may be indicative of individuals with OA in older populations, particularly males. Our metabolite signature modeling method is likely to increase classification power in validation cohorts.


Introduction
Commencing at age 50, there is a steep increase in the incidence of symptomatic osteoarthritis (OA) and the number of individuals undergoing total knee replacement (TKR) [1][2][3]. Individuals of high body mass index (BMI) and female sex also have increased risk of TKR due to primary OA [4,5]. The World Health Organization defines individuals with BMI ! 30 as obese. There are currently no effective biomarkers to identify individuals with advanced OA.
The metabolome represents the cumulative output of metabolic processes occurring within an individual and includes compounds such as lipids and amino acids, among others. In addition to the amino acid arginine [6], select lysophosphatidylcholine (lysoPC) to phosphatidylcholine (PC) analogue ratios are altered in plasma from patients with OA compared to healthy adult volunteers (HV) and the ratio of total lysoPCs to PCs are predictive of TKR in 10 years follow-up [7]. It is unclear, however if a signature of individual metabolite levels, specifically of lysoPC or PC analogues, in a combined signature, is capable of classifying OA from healthy individuals. In addition, optimal methods to help improve predictive metabolite selection and prediction modeling in cross-sectional cohorts to aid in successful prediction in external validation cohorts have not been investigated.
In this study, we sought to determine if a signature of metabolite levels could be predictive of OA vs healthy volunteers using a selection and modeling method to improve selection and predictive consistency. Using plasma from a cohort of patients undergoing TKR surgery due to primary OA and HV, we measured lysoPC and PC analogues. Stratifying along age, sex and BMI, we determined unique signatures based on each stratum using systematic univariate modeling followed by prediction analysis using various multivariate modeling algorithms for comparison. Thus, we present a method to investigate demographically-stratified populations from a single cross-sectional cohort to determine the best possible predictive metabolite signature of individual metabolite levels that can be used in future validation studies.

Study participants
Patients receiving TKR due to primary OA were obtained from the Newfoundland Osteoarthritis Study (NFOAS) [8]. HV were from the The Complex Diseases in the Newfoundland population: Environment and Genetics (CODING) study [9]. Both OA patients and HV were from Newfoundland & Labrador, Canada. Knee OA diagnosis was made based on American College of Rheumatology clinical criteria for classification of idiopathic OA of the knee [10] and judgment of attending orthopedic surgeons. Controls were those without an OA diagnosis in any joints based on medical information collected by a self-administered questionnaire. The distribution of clinical and demographic variables in strata and a comparison between OA and HV individuals is shown in Table 1.

Metabolite profiling, signature determination, and predictive modeling
Blood samples were collected after minimum 8 hours of fasting. Blood was collected into K 2 EDTA-plasma tubes. EDTA-plasma was separated from whole blood by centrifugation at 1500 rcf for 10 mins at 4˚C. Plasma was aliquoted at stored at -80˚C until use. Plasma was thawed on ice and metabolite profiling was performed using Liquid Chromatography(LC)/ Mass Spectrometry (MS)/MS using the Waters XEVO TQ MS Ultra Performance LC/MS/MS system (Waters Limited, Mississauga, Ontario, Canada) coupled with Biocrates AbsoluteIDQ p180 kit. All analytical metabolite quantification was performed using the Absolute IDQ- coupled MetIDQ software package (Biocrates Life Sciences AG, Austria), as described [8].
Only lysoPC and PC analogue concentrations were used for analysis in this study. Metabolite concentrations were adjusted for batch effects and log (x+1) transformed to normalize variable distributions for modeling using linear effects. The cohort was stratified based on age ( 50 or > 50 years), BMI (! 30 or < 30 kg/m 2 ), and/or sex. The information for each stratum is presented in Table 1.
Subpopulations (OA vs. HV, stratified by age) were assessed for metabolite variance homogeneity to check for structural differences. Homogeneity was tested using PERMDISP2 [11] and significance was measured using Tukey's HSD test. For model building, we pre-selected individual metabolites within each of the stratified populations using predictive area under the curve (AUC) in a logistic regression using a non-parametric bootstrap [12]. We randomly sampled N individuals with replacement to generate 1000 training sets to which we fit a logistic regression with single metabolites, and then estimated predictive AUCs on individuals not included in the corresponding training set. Under the bootstrap, this is about 1/3 of the samples, constituting the test set. Using the replicates, we generated an empirical distribution of AUCs for each metabolite and selected those with consistent predictive power (AUC > 0.5 at the 2.5% quantile).
Metabolites above the cutoff for each stratified population were used as inputs into three predictive algorithms for classification: principal component regression with logistic regression (PCR), partial least squares regression (PLS) with logistic regression, and simple logistic regression. PLS and PCR approaches were used for dimensionality reduction as metabolites tended to be correlated. The use of multiple metabolites projected onto components of correlation for PLS and variation for PCR is to achieve a more robust signal, compared to reducing metabolite lists to achieve a higher AUC, which may be specific to the dataset and thus overfit the model. For each modeling algorithm, the bootstrap process was repeated 1000 times to generate an empirical distribution of AUCs. For PCR and PLS, the number of components was selected by minimizing overfit, as measured by the differences between average training and test set AUCs from the bootstrap. The results indicated that the first principal component was optimal in both models. Aggregate concentrations of metabolite groups [(lysoPC, diacyl PC (PCaa), acyl-alkyl PC (PCae)] from all measured metabolites were also used in separate logistic models. The process was repeated for all strata. The modeling process is presented in Fig 1. We have also provided the statistical code used for analysis using R package (https:// www.r-project.org; S1 Document).

Study Approval
The study was approved by the Health Research Ethics Authority (HREA) of Newfoundland & Labrador (reference number 11.311) and written informed consent was obtained from all participants prior to inclusion in the study.

Results
The cohort utilized for the study is presented in Table 1. A total of 346 individuals (152 OA, 194 HV; 159 males, 187 females) were studied. The OA population was older and had higher BMI (all P < 0.05), but the groups had similar sex distributions. Since age was significantly different, our study focused on individuals over 50 years old. This reduced the mean difference in age between OA and HV to 4.5 years compared to 13.6 years in the entire cohort (Table 1). Furthermore, variance of metabolite levels between age-stratified cohorts was significantly different (P < 0.05) whereas the variance of metabolite levels between OA and HV individuals over age 50 was not significantly different (P > 0.05), further suggesting a need for age stratification. The variance of metabolite levels between OA and HV individuals in BMI-stratified groups, either in the entire cohort or in individuals over the age of 50, was not significantly different (P > 0.05). This suggested that age, rather than BMI was the major confounder, even though BMI between HV and OA individuals was significantly different in stratification of age > 50 years with BMI ! 30 or BMI < 30 (Table 1).
Qualitatively, we were unable to identify overt trends between males, females, HV and OA patients over 50 years old using heat-map analysis (Fig 2A), suggesting that smaller changes in individual metabolite levels are important for detecting OA. We therefore analyzed each stratum population to identify individual metabolites predictive of OA using univariate receiver operator AUC analysis (Table 2). A list of predictive metabolites for OA was identified using bootstrapped logistic regression analysis (Fig 1). Most of the strata generated unique signatures of OA-indicative metabolites, except for females and individuals with BMI < 30 kg/m 2 , which did not have any metabolites that met our selection criteria ( Table 2).
Comparing signatures obtained from stratum age > 50 years old to the same aged populations also stratified based on sex, we determined that the signature identified in males age > 50 years old encompassed all of the metabolites identified within the age > 50 years population, along with the single metabolite identified for the age > 50 years female population ( Fig 2B, Table 2). The signature from males age > 50 years old stratum also encompassed all but 2 metabolites derived from the age > 50 years with BMI ! 30 or < 30 kg/m 2 stratified populations ( Fig 2C). This suggests that lysoPC and PC analogues have most power to detect OA in older males.
Next, we investigated multivariate algorithms to identify ones with the most consistent classification power (Fig 1, Table 3). Analyzing differences in AUC we found that multivariate logistic modeling was poor at maintaining predictive consistency. For principal component with logistic regression (PCR) and partial least squares with logistic regression (PLS) approaches on the top component, PCR was slightly more consistent than PLS modeling based on absolute mean differences, being more consistent at the 2.5% and 50% quantiles, but slightly worse at the 97.5% quantile. PCR modeling also tended to be more conservative, as training set AUCs were marginally smaller compared to PLS (Table 3). Based on PCR modeling, we found that the males > 50 years old signature was most accurate in differentiating OA and HV, with a median AUC of 0.751 and 0.752, for the training and test sets respectively ( Fig  2D, Table 3).
Finally, we evaluated whether all measured phospholipids aggregated by type, namely lysoPCs, PCaas, and PCaes, could be used as predictors of OA ( Fig 2E, Table 3). Consistent classification accuracy above random was not achieved for all strata using the same modeling algorithms at a 2.5% quantile cutoff. For PCR and PLS modeling, equivalent quantiles of AUC were between 5.2-20.8% and 4.9-18.4% lower, respectively, compared to using signature metabolites. Thus, performance was inferior to using signature metabolites modeled using PLS or PCR (Table 3).

Discussion
We have outlined an approach to identify individual metabolites strongly indicative of OA and a method of using them in classification modeling. Our results showed OA vs HV classification using metabolite levels was strongest in older males, resulting in a signature of 32 metabolites that individually had power to classify OA vs HV. Furthermore, we determined that modeling the signature using PCR with a single component resulted in the most consistent classification accuracy, as measured by mean AUC differences between bootstrap training and test sets. Finally, we concluded that using a signature-based list of metabolites and PCR/PLS modeling https://doi.org/10.1371/journal.pone.0199618.g002 Table 2. Metabolites with a 2.5% quantile area under the receiver-operator curve ! 0.5 determined by bootstrapped logistic regression of the stratified study population described in Table 1   AUC values were generated using partial least squares (pls), principal component analysis and logistic regression (pcr) and multivariate logistic regression (log) alone. Aggregated lysphosphatidylcholine, diacyl-phosphatidylcholine (PCaa) and acyl-alkylphosphatidylcholine (PCae) concentrations were also modelled (sum) in the same manner for each stratified group.
Differences in training and test set were calculated and the mean absolute difference across all stratified groups was calculated for each quantile and for all quantiles to identify the model with least amount of overfitting between bootstrapped test and training sets. Age (in years), body mass index (BMI; in kg/m 2 was more effective than aggregate modeling of all metabolites with the same methods, suggesting metabolite subsets are superior for detecting OA. The statistical approach for this study (Fig 1) was chosen to reduce overfitting, allowing for greater confidence in the signatures we have identified and their potential classification power in independent datasets. With limited metabolomics research conducted on circulating serum or plasma from healthy and OA cohorts to date, validation has yet to be conducted on external cohorts. This is the first method in OA metabolomics research, to our knowledge, to focus on identifying cohort-specific signatures via demographic stratification. Independent cohort data will be necessary to conclusively determine if the modeling algorithm and identified signatures are universally predictive of OA.
Since metabolite signatures were strongest in the strata of males age > 50 years, metabolic processes may be more consistently modulated between OA and non-OA individuals in this stratum. Thus, in males, lysoPCs and PCs may be indicative of a response to OA pathophysiology or a direct contribution to symptoms/pathogenesis. The lack of this signal in females, given their higher general population incidence, may suggest varying etiology, a more biochemically heterogeneous disease population where metabolite signatures are less pronounced. However, in all of our metabolite signatures, confounding clinically-relevant variables were not controlled for in this study due to the lack of information from HV, a limitation of this study which may affect both signature elements and predictive outcomes. Controlling these confounding clinical variables using the method described herein could also increase both accuracy and precision of OA prediction.
We identified that compared to aggregates of metabolite types [7], select metabolite signatures are capable of classifying OA vs. HV in select strata. We also determined that this was as good as or better than aggregate levels alone, suggesting that select metabolites likely drive classification of individuals with OA. This is consistent with studies in other diseases where metabolite signatures were capable of classifying individuals with various pathologies including high-altitude pulmonary edema [13], multiple sclerosis [14] and pediatric tuberculosis [15].
Our study does however, have some limitations based on the clinical data that was available for use. For instance, we could not stratify based on metabolic disorders, including hyperlipidemia due to the lack of available clinical data, which could alter metabolite levels in the blood and be associated with hyperlipidemia as opposed to OA directly. Hyperlipidemia is a risk factor for OA [16] whereas metabolic syndromes are correlated with knee OA prevalence in select populations [17,18]. Furthermore, although the incidence of diabetes mellitus is known in our cohort, we did not stratify or exclude individuals based on this factor, as we were unaware if this was well controlled by medication or lifestyle. Removing these individuals would also reduce the power for some of our strata. However, we found 4 metabolites that were significantly different in their levels between individuals with or without diabetes mellitus (DM) within the age > 50 years (3 metabolites) and males age > 50 years (1 metabolite) strata that overlapped with metabolites identified in our signatures from Table 2 (S1 Table). Removing these metabolites from our signatures and running PCR modeling of the remaining metabolites in the two affected strata resulted in minimal changes to the resulting AUCs (S2 Table). Curiously, there were more metabolites significantly different in the HV compared to the OA populations when comparing individuals with or without DM, suggesting that individuals with OA are more similar, even if they have DM (S1 Table). Overall, this suggests that the presence of DM is not a major confounding variable in our study.
Thus, we determined that individual lysoPC and PC analogues in plasma were collectively able to detect OA with some accuracy. Inclusion of other circulating biomarkers, like amino acids [6], cytokines [19] and microRNAs [20], into modeling algorithms may improve prediction.
Supporting information S1 Document. R v3.1.1 coding used for dispersion tests, univariate AUC analysis and, multivariate regression AUC analysis. (DOCX) S1 Table. Metabolites significantly different between individuals with or without diabetes mellitus within healthy control volunteers or individuals with osteoarthritis stratified by age, body mass index (BMI; kg/m 2 ) and sex. (DOCX) S2 Table. Quantiles of the principal component regression with logistic regression area under the multivariate receiver operator curve values of signatures from Table 2 with significantly altered metabolites identified in diabetic vs non-diabetic individuals (S1 Table) removed from the signatures. (DOCX)