A weighted relative difference accumulation algorithm for dynamic metabolomics data: long-term elevated bile acids are risk factors for hepatocellular carcinoma

Dynamic metabolomics studies can provide a systematic view of the metabolic trajectory during disease development and drug treatment and reveal the nature of biological processes at metabolic level. To extract important information in a systematic time dimension rather than at isolated time points, a weighted method based on the means and variations along the time points was proposed and first applied to previously published rat model data. The method was subsequently extended and applied to prospective metabolomics data analysis of hepatocellular carcinoma (HCC). Permutation was employed for noise filtering and false discovery rate (FDR) was used for parameter optimization during the feature selection. Long-term elevated serum bile acids were identified as risk factors for HCC development.

behavior and curve-fitting to analyze the compounds. Although these methods are compatible with short time-series datasets, each observed time series is assumed as a smooth random curve. However, when dealing with detailed time-series data where specific time points must be treated differently, corresponding data processing methods are needed.
Hepatocellular carcinoma (HCC) is one of the most lethal cancers 28,29 , and its incidence and mortality rates continue to increase 30 . However, the mechanism of hepatocarcinogenesis remains obscure because of the complicated interactions of multiple factors and individual genetic variations, impeding early clinical intervention before the development of HCC. Relatively effective treatments are available when HCC is diagnosed early. HCC patients often have a history of chronic liver diseases, leading to the introduction of screening programs among high-risk populations 31 , such as those infected with hepatitis virus B (HBV) in Qidong, China (a high-incidence area of HCC due to the high prevalence of HBV infection), who undergo HCC screening every half year. In addition, a sample library has been established in Qidong for HCC pathogenesis and early diagnosis studies [32][33][34] .
In this study, a weighted relative difference accumulation algorithm (wRDA) and its extended form were proposed. The wRDA method was first used to treat our previously published rat model data, and its extended form was further applied to a prospective cohort study of HCC patients with the aim of revealing earlier HCC diagnosis biomarkers and metabolic dysregulations contributing to hepatocarcinogenesis.

Results
The application of the wRDA to metabolomics data from the rat HCC model. The proposed wRDA was first applied to our previously published data for a rat HCC model induced by diethylnitrosamine (DEN) administration 35 . In that study 35 , 52 differential metabolites were identified, of which three, taurocholic acid (TCA), lysophosphoethanolamine 16:0 (LPE 16:0) and lysophosphatidylcholine 22:5 (LPC 22:5), were defined as ''marker metabolites'' for distinguishing the different stages of chemical hepatocarcinogenesis. LPE 16:0 and TCA were more discriminative between the disease group and control group, whereas LPC 22:5 was more discriminative between the HCC and non-HCC samples.
Parallel to our previous feature-defining process, a two-level data analysis procedure employing the wRDA (Fig. 1) was performed to select meaningful features to discriminate between the models and control, and between HCC and non-HCC samples. In the first level, 152 ion features were removed by means of permutation, leaving 1092 features with a false discovery rate (FDR) of 0 to constitute feature subset 1. Then, Support Vector Machine 36 (SVM) was conducted based on the top 20 features ( Fig. 2A) ranked by the wRDA. Five fold cross validation was conducted 50 times. The average accuracy rate was 98.85% 6 0.66%, demonstrating that the top ranked variables have a strong ability to distinguish disease samples from controls. These metabolic features include two bile acids (TCA and tauroursodeoxycholic acid (TUDCA)), LPCs, and LPEs with different acyl chains and unsaturation levels. These results indicate a disturbance of lipid metabolism in DEN-induced liver disease.
In the second level, the features in feature subset 1 were analyzed again to calculate their ability to characterize the metabolic status of the three different liver diseases. Using the top 20 ranked variables (Fig. 2B), the average accuracy rate of the SVM classifier for discriminating HCC and non-HCC samples was 93.12% 6 1.92%, demonstrating that informative features can be defined by weighting the time points according to their importance (here the importance of the time points was decided according to prior knowledge). The combinations of these features can denote disease progression towards HCC.
To investigate the discrimination abilities of the top ranked 20 features defined by wRDA to distinguish the liver diseases from the controls, their receiver operator characteristic curves are further drawn with their area under curves (AUCs) calculated. Among 15 features (after deleting the redundant ions from the same metabolite), three features TUDCA, feature with m/z of 636.3415 and LPE   35 , here further explanation is omitted. Collectively, the metabolic features derived from the rat model metabolomics data demonstrate the excellent performance of the wRDA in feature-ranking and demonstrate its potential in analyzing time-series metabolomics data.
The application of w 2 RDA in HCC prospective metabolomics data analysis. In this prospective cohort study, samples of 5 time points from 11 HCC cases were collected during screening over 3.5 years. The study aimed to identify prospective features of HCC, but the biological events in the other time points before HCC occurrence are unknown. Moreover, in contrast to animal model samples, the collection times of samples from patients in the same stage were not uniform. Therefore, parameter k was introduced into the wRDA to reduce the influence of sample storage time. The settings of v i for each time point and k j for each sampling time affect the measurement of the features. To define the most suitable settings for v and k, linear function, exponential function and proportional function were tested with different changing factors (equal weights were included as a special case of linear function). For n top-ranked features, the lowest FDR was adopted to evaluate the settings of k and v. As shown in Table 1, when v was an exponential function, the minimum values of the lowest FDR were obtained for each k's function setting. When k was also an exponential function, the minimum FDR values in each column were derived for the vast majority of different function settings of v for n decreasing from 50 to 30. Furthermore, 194 k-v pairs were derived from different functions and their corresponding changing factors, resulting in a low FDR of less than 5% (including 192 k-v pairs with the lowest FDR equaling 0 at n 5 30 and 2 k-v pairs with the lowest FDR equaling 2.86% at n 5 35). Among the 194 k-v pairs, the probability of both v and k being exponential was the highest (25.26%). Therefore, exponential function is more suitable for k and v than other functions.
When k and v were fixed as exponential functions, the lowest FDR was 0% at n 5 30, and there were 47 paired values of changing factors for k and v (Table S1). To minimize the weight differences among the sampling points, the smallest changing factor of k was chosen, q k 5 0.5. Once k was selected, the smallest changing factor of v was defined as q v 5 0.6.
Under the optimized settings for k and v, the features were ranked according to their w 2 RDA scores. The top 30 ranked variables ( Table 2) were chosen as the most important metabolic features related to HCC development.
The most important variables were bile acids, with the exception of dihydroxyandrostenone sulfate and LPE 18:2. The serum concentrations of the primary bile acids cholic acid (CA) and chenodeoxy-   cholic acid (CDCA) were significantly higher in the HCC group than in the control group (at T 0 ) when malignant hepatic tumors were identified. Four other bile acids, the secondary bile acids deoxycholic acid (DCA) and taurodeoxycholic acid (TDCA) and the sulfated bile acids glycodeoxycholate sulfate (GDCAS) and glycochenodeoxycholate sulfate (GCDCAS), were also elevated in HCC sera compared to controls. When the serum levels of bile acids were compared over the entire monitoring time period, most were significantly elevated in patients in which HCC occurred compared to individuals that were hepatitis B surface antigen positive (HBsAg 1 ), with the exception of hyodeoxycholic acid (HDCA) ( Table 2). The relative serum levels of bile acids at each time point are shown in Fig. 3.  Discussion HCC usually develops from chronic liver diseases, and a long time is required for the formation of malignant hepatic tumors. The mechanism underlying the occurrence and development of HCC remains to be elucidated. Bile acids are synthesized in the liver and their functions are not limited to facilitating the absorption of lipids and lipid-soluble nutrients 37,38 but also include acting as signaling molecules to regulate glucose and lipid metabolism 39,40 and apoptosis 41 . Bile acids are detergents and are cytotoxic, and their concentrations are tightly regulated under normal physiological conditions 42,43 . We previously demonstrated that glycocholic acid (GCA) and glycochenodeoxycholic acid (GCDCA) are elevated in hepatitis, cirrhosis and HCC accompanied by cirrhosis 44 . The serum levels of seven bile acids were quantitatively measured and compared among HCC patients without liver cirrhosis and hepatitis, HCC patients with liver cirrhosis and hepatitis, benign liver tumor patients with liver cirrhosis and hepatitis, and healthy controls, and elevated serum GCDCA, GCA and TCA levels and decreased serum CDCA levels were correlated with liver cirrhosis and hepatitis 45 . We previously demonstrated that conjugated GCA, GCDCA, TCA and taurochenodeoxycholic acid (TCDCA) are potential biomarkers of liver cirrhosis 46 .
Fasting serum levels of primary bile acids 47 can be affected by enterohepatic circulation 48 , leading to intra-individual variations. Thus, multiple time points of circulating bile acids were compared together within the monitoring period instead of at a single time point. This comparison revealed that all bile acids except HDCA were significantly higher in HCC patients than in HBsAg 1 controls. The more hydrophobic secondary bile acids DCA and lithocholic acid (LCA) have been reported to increase HCC risk 49,50 .
Activation of the YAP pathway was recently shown to be responsible for bile acid-dependent tumor promotion 51 . The development of spontaneous liver tumors in a Fxr -/-Shp -/double-knockout (DKO) mouse model was employed in the study to produce chronically elevated bile acid levels, which enabled the study of the mechanism of hepatic malignant tumor promotion by long-term high circulating levels of the bile acids CA and CDCA in mice 51 . Compared to the HBsAg 1 control group, serum CA, CDCA and DCA levels were slightly elevated in the HCC group during the two-year monitoring period, and their glycine-and taurine-conjugated forms were elevated to a greater extent. Thus, the increased serum levels of bile acids may be due to leakage from damaged hepatic cells or the alteration of bile acid transfer protein activity 52 rather than upregulation of bile acid synthesis. When more bile acids enter the blood, they may further intensify hepatic injury because of their cytotoxic nature and may simultaneously act as signaling molecules to promote hepatic tumor formation.
No differences in levels of ursodeoxycholic acid (UDCA), which has been reported to have protective effects 53,54 , were observed in HCC patients and HBsAg 1 controls during the two years before or at HCC diagnosis, whereas levels of its taurine-conjugated form were slightly elevated in HCC patients. The sulfated bile acids GDCAS and GCDCAS were elevated in HCC patients during the two years prior to diagnosis. Sulfotransferase-2A1, which has been reported to be underexpressed in HCC tumor cells 55 , catalyzes the sulfation of bile acids for their elimination and detoxification 56 . Sulfotransferase activities have been reported to decrease with the severity of liver disease from steatosis to cirrhosis 57 . The long-term increase in sulfated bile acids in HCC patients may be due to their increased availability for sulfation rather than enhanced SULT2A1 activity.
Because chronic hepatitis and cirrhosis are typically precursors of HCC, in combination with the above evidence that bile acids promote hepatic tumor formation, it is reasonable to speculate that long-term high circulating bile acids are potential high-risk factors for HCC. TCA is elevated since week 6 in model rats treated with DEN compared to controls 35 . The 13 bile acids mentioned above were extracted from the DEN-induced rat HCC model data acquired in positive mode, which revealed that the 10 bile acids (CA, CDCA, DCA, GCA, GCDCA, glycodeoxycholic acid (GDCA), TCA, TCDCA, TDCA and tauroursodeoxycholic acid (TUDCA)) detected were elevated in sera in the model group compared to the control group since week 6 ( Supplementary Fig. S2), coincident with the appearance of hepatic cell injury due to DEN treatment. It has been reported that 40% of HCC patients infected with HBV from Qidong have high exposures to aflatoxin B1 58 . Although the mechanism of HCC pathogenesis may vary greatly because of different etiological agents, the common long-term elevated serum bile acids were observed before HCC occurrence. Collectively, it can be speculated that the elevated levels of circulating bile acids in chronic liver disease may play an important role in the process of malignant hepatic tumor formation in both humans infected with HBV and DEN-treated rats.
In this article, to identify discriminative metabolites that may reflect dynamic biochemical developments, a wRDA method based on the weighted mean and variance analysis along the time points was proposed. Rather than screening differentially expressed variables at isolated time points as in static methods, the wRDA can investigate variables comprehensively along all time points in feature selection. Moreover, weighting the time points emphasizes the influence of relevant, important time points by setting relatively larger corresponding weights, and vice versa. Thus, high efficiency can be achieved in identifying the key differences between two groups along the entire time course.
The application of the wRDA to the rat metabolomics data demonstrated that it is an effective method for defining metabolic features that may be related to disease status. In the prospective cohort study, continuously high serum bile acid levels within a two-year monitoring time period were identified as risk factors for HCC development. The weights of the time points can be decided based on prior knowledge or optimized by the lowest FDR. Other functions and changing factors for weighting can be simulated, and other optimization standards can be introduced in further studies. Collectively, our proposed method, by analyzing the weighted relative difference accumulation along the time dimension, effectively defines the features of dynamic metabolism related to disease development.
Methods wRDA and its application to the rat liver disease model. wRDA. When analyzing the metabolomics time course data of two different groups, for simplicity, let C denote the control group and M denote the model group. Let T i denote a time point, 0 # i , N, where N is the number of the time points.
In bioinformatic data analysis, the method used to measure the discriminative ability of a feature among different groups is a key consideration. SAM 59 scores the ''relative difference'' of a gene over repeated measurements according to the mean and the standard deviation. Based on the idea of ''relative difference'', the wRDA considers the variations of the means along the time points to dynamically study the biological process and screen biomarkers that reflect differences between the two groups and characterize the development of the model group. A variable with a higher wRDA score is more discriminative between the two groups. The detailed principles of the wRDA are outlined as follows: In metabolomics time-series studies, metabolites are measured at each time point. Because differences in metabolite levels in the two groups may occur in a process, the accumulation of distance between the mean values of a given feature f at all time points, D(f), reflects the discriminative ability of feature f: where m C,f (i) and m M,f (i) are the mean values of feature f at time point i in the C and M groups, respectively and v i represents the weight of time point i. Different time points may play different roles in the development of differences between the two groups, resulting in different weight settings. In particular, some time points may occur at typical stages for biological events and hence merit greater attention and relatively large weights. Furthermore, the standard deviation is applied to enable a fair comparison among the discriminative abilities of features. Let www.nature.com/scientificreports SCIENTIFIC REPORTS | 5 : 8984 | DOI: 10.1038/srep08984 where s C,f (i) and s M,f (i) are the standard deviations of feature f at time point i in the C and M groups, respectively. Hence the ''weighted relative difference accumulation'' of a feature f between the C and M groups over all time points is calculated as wRDA(f): : ð3Þ e is a small positive value introduced to moderate or regularize the wRDA score and potentially reduces the relative impact of small variances. In this study, e was set to 0.005.
Data source for the rat liver disease model. Metabolomics dynamic data for the rat liver disease model 35 were employed using the model obtained from the Shanghai Experimental Animal Centre. Serum samples were collected from two groups, control rats and rats with liver disease induced by DEN, every two weeks from week 6 to week 20. A total of 8 monitoring time points were obtained. In this animal model, the serial progression of hepatocarcinogenesis includes three disease stages: the inflammation stage (week 6-week 8), the cirrhosis stage (week 12-week 14), and the HCC stage (week 18-week 20). All samples at each time point were collected synchronously in this animal experiment.
The application of the wRDA to the rat liver disease model. The wRDA was first applied to the metabolomics data of the rat liver disease model. First, features whose values equaled zero in more than 20% of the samples 60 at each time point were removed, leaving 1289 features. Then, outlier correction was conducted (a sample for feature f is an outlier if its value is beyond the range m(f) 6 2s(f), where m(f) is the mean value of f and s(f) is the standard deviation of f in the corresponding group 13 ). Assuming that the feature followed a prior distribution, the outlier was replaced by a random selected sample value following this distribution.
To select the features reflecting the different developments between the two groups and define the features discriminating HCC samples from non-HCC samples, the wRDA was applied at two levels ( Fig. 1).
(1) In the first level, to study the dynamic differences between the liver disease group and the control group, equal values of 1/8 were assigned to v at all time points. The features defined together reflect the entire liver disease status of the model group from week 6 to week 20. Permutation was conducted 200 times to filter noise and noninformative features.
(2) At the second level, the wRDA was applied again to measure the variables in feature subset 1 (Fig. 1), and the weights (v) of three time points, week 6, week 12 and week 18, which were defined as the typical stages of hepatitis, cirrhosis and liver cancer, were set to 0.3, 0.3, and 0.4, respectively. The stage at which liver malignant tumors occur is the most important and merits greater attention, as reflected by a larger weight. The weights (v) of the other time points were set to zero. The weight adjustment allowed the most informative features characterizing the three different liver diseases to be targeted, particularly those capable of distinguishing HCC from non-HCC. w 2 RDA and its application in the HCC prospective cohort study. w 2 RDA. In epidemiological screening, people with high risk are checked at certain time intervals. One time point may lie in the onset stage of a disease or a malignant tumor, and other time points may be a certain time interval before or after the key biological event. However, not all patients with a disease or a malignant tumor were spot at a uniform screening period, and thus samples at each disease stage may be collected at different sampling times. To measure the features more accurately, a weight k for different sampling times was introduced: where p i is the number of the sampling times at time point i; m C,f (i, j) and m M,f (i, j) are the average levels of feature f at the sampling time j of the i th time point in the C and M groups, respectively; and k j is the weight of the j th sampling time. The extended standard deviation by considering the sampling time differences at each time point is defined as follows: where s C,f (i.j) and s M,f (i,j) are the standard deviations of feature f at the sampling time j of the i th time point in the C and M groups, respectively.
HCC prospective cohort study and data collection. Serum specimens were obtained from a prospective cohort in Qidong, Jiangsu Province, China. From May 2009 to October 2012, residents of the Qidong area were invited for a health examination as part of the HCC screening study. Table 3 shows the baseline characteristics of the enrolled HCC and HbsAg 1 control subjects. Each participant underwent serological hepatitis tests (HBsAg, hepatitis B e antigen (HBeAg), anti-hepatitis C virus (HCV)), abdominal ultrasonography (US), serum alpha-fetoprotein (AFP) and alanine aminotransferase (ALT) tests at approximately six-month intervals. If an individual had abnormal results from US or higher AFP levels (greater than 20 ng/mL), then intensive surveillance by computed tomography (CT), magnetic resonance imaging (MRI), and/or hepatic angiography were employed to identify the space-occupying lesion. In total, 11 HCC cases were identified in a 3.5-year screening period. Their serum samples at the time of HCC diagnosis and samples within the 2 years preceding the initial HCC diagnosis were selected for analysis. Another 22 HBsAg 1 individuals were taken as a positive controls with matched age, sex, sample collection time points and storage conditions (220uC). The details of serum sample preparation, liquid chromatography-mass spectrometry (LC-MS)-based metabolic profiling and data preprocessing are available in the supplementary material.
The application of the w 2 RDA in the HCC prospective cohort study. For simplicity, the HCC and HBsAg 1 control groups are also denoted as M and C, respectively. T 0j , T 1j , T 2j , T 3j and T 4j (Fig. 4) are used to mark each stage (time points). T 0j represents the stage when HCC was diagnosed, whereas T 1j , T 2j , T 3j and T 4j represent the stages 0.5, 1, 1.5, and 2 years before T 0j , respectively. HCC patients were identified from the participants at four screening periods in this study; thus there are four sampling times for each stage (Fig. 4). Taking the T 0j stage as an example, patients were diagnosed with liver cancer in May 2011, Nov. 2011, May 2012 and Oct. 2012. The corresponding serum samples at 6 months, 12 months, 18 months and 24 months prior to HCC diagnosis were then collected. The longest time interval of the samples at each stage from different sampling times was 1.5 years. Therefore, the w 2 RDA was applied to further consider the possible influence of different sampling times. The aim of the Qidong cohort study was to identify prospective features related to the occurrence of HCC. The settings of the weights v and k could affect the feature measurement of w 2 RDA. T 0j is the onset stage of HCC and therefore is the most important. The closer the other time points are to T 0j , the more similar their metabolic characteristics are to those of HCC. Hence, for v (and k), three different setting methods were tested: a linear function, a proportional function and an exponential function. FDR 59,61 was adopted to evaluate the results under different weight settings. For the linear function v i 5 1.0 1 (N-i-1)q, 0 # i , N, q is a parameter factor. When q 5 0, all time points have the same weight. For the proportional function v i 5 (1.0 1 q) (N-i-1) , 0 # i , N. Many metabolomics experiments use natural exponential functions to estimate time-varying profiles 6 ; in this case, v i 5 e (N-i21)q , 0 # i , N. In the latter two functions, q is also a changing factor. In the Qidong cohort study, to consider the influences of all the monitoring time points on metabolism, the differences of v i should not be too large. q was restricted from 0.1 to 1.0 and was tested with a step increment of 0.1. In addition, T 0j is the most important stage when malignant hepatic tumors are discovered, and thus its weight should be larger than those of the others. Similarly, the weight k was also tested using a linear function, proportional function and exponential function. The sample with the shortest storage time should have the greatest weight.