Hybrid feature detection and information accumulation using high-resolution LC- MS metabolomics data

The overall workflow

The workflow is illustrated in Figures 2. The hybrid peak detection approach involves the building of a known-feature database. The information can come from two sources – databases of known metabolites and historical data. The historical data must be generated from the same LC/MS procedure as the data being analyzed.

Once the database is available, our method employs a four-step procedure to detect features from a new batch of LC-MS profiles. First, it analyzes the data using traditional unsupervised approach with relatively stringent criteria for peak calling ⁹. A feature list is generated. Second, the newly found features are compared with those in the database. If a feature found in the new data matches a known feature closely, it is identified as the known feature. Newly found features that do not match to any known feature are considered as potential new features. Third, we combine all the features in the database and the potential new features into a candidate feature list. For every candidate feature in this list, we conduct targeted peak detection to search the LC-MS profiles. In this step, a peak filter that is less stringent than in the first step is employed, as we have higher confidence of a peak’s existence. Fourth, after targeted peak detection, the database is updated to include new information, and a feature table is generated for the current data for downstream analysis.

Initializing the database

The apLCMS package utilizes the database in R data frame format. Key components of the database include feature ID, possible chemical identity, m/z (theoretical value for known metabolites or mean observed value for unidentified historically observed features, standard deviation, range), retention time (mean, standard deviation, range), log intensity (mean, standard deviation, range), and the frequency the feature is detected in historical data (Supporting Figure 1). Features derived from known metabolites may not have retention time and intensity information available, while features from historical data may not have their chemical identities available.

Various methods can be used to store the database. Simple ones include tab-delimited text files and directly saving the data table in R binary format, both of which can be done in a single line of command. We illustrate them on the apLCMS website. More complex methods include using the R interface to relational databases, such as the RMySQL package available from CRAN. The database itself may store more information than apLCMS utilizes, and only necessary information needs to be extracted for use in apLCMS.

In the current study, we merged metabolites from HMDB ²⁹ and features found from the example dataset of the apLCMS package ^{9, 20}. The data was generated using anion exchange column with formic acid gradient combined electrospray ionization (ESI) ³⁰. Briefly, human plasma samples were collected, and protein precipitation was performed by adding acetonitrile (2:1, v/v). Analyte was separated using the Hamilton PRPX-110S (2.1×10 cm) anion exchange column together with a short, end-capped C18 pre-column (Higgins Analytical Targa guard) for desalting and optimal separation. A formic acid gradient was used. Solution A was made of 0.1% (v/v) formic acid in a 1:1 water:acetonitrile mix. Solution B was made of 1% (v/v) formic acid in a 1:1 water:acetonitrile mix. Mass spectrometry data was collected using either a Thermo LTQ-FT mass spectrometer or a Thermo Orbitrap-Velos mass spectrometer (Thermo Fisher, San Diego, CA) using and m/z range of 85 to 850. For more details please refer to Johnson et al ³⁰.

Some metabolites have the same chemical composition, hence identical m/z values. For a group of metabolites sharing m/z, a single entry was created, which contains several possible chemical identities. When new information becomes available to separate those metabolites, the single entry can be split into multiple entries. To unravel such groups needs techniques such as LC-MS/MS, which is out of the scope of this study.

At the 10 ppm tolerance level, we merged known metabolites with features found in the example dataset. For HMDB metabolites, we only considered one ion form [M + H]⁺, because ESI mostly generate protonated molecular ions. The choice of ion form is highly dependent on the experimental platform. Users can incorporate other ion forms, as well as isotopic forms of the known metabolites by generating entries based on the neutral mass of the metabolites and creating their own database. Detailed format of the database is explained in the manual on the apLCMS website. Although other ion forms are ignored, they can be incorporated into the database if they are found in real data. From features found in the example dataset, we only incorporated reliable ones into the database by requiring a feature to be present in at least 50% of the profiles.

Comparing identified features in new data with the database and generating candidate feature list

Once a feature list is generated from the new dataset using traditional unsupervised peak detection and peak alignment, we need to merge it with the features in the database by finding the correspondence between the newly detected features and those in the database. The apLCMS package conducts retention time deconvolution using the bi-Gaussian mixture model to separate features sharing m/z ²⁰. Thus the feature list from the new data may contain features with almost identical m/z values and different retention time values.

First, we match the two lists at a certain m/z tolerance level ϕ_m, e.g. 10 ppm, and retention time tolerance level ϕ_r which can be either user input or found by adaptive search from the data ⁹. Second, in the case of multi-to-one or multi-to-multi match, we take the features involved, and find the distance matrix D with newly found features in the rows and known features in the columns. Every element of D is computed using the Mahalanobis distance with $\left(\begin{array}{cc}{\varphi }_m^2& 0\\ 0& {\varphi }_r^2\end{array}\right)$ as the variance-covariance matrix. Then we iterate the following procedure: (a) find the smallest element of D and its location in the matrix; (b) declare the pairing between the new feature of the corresponding row and the known feature in the corresponding column; (c) remove the corresponding row and column; (d) repeat the above three steps until the smallest element is larger than a pre-specified threshold, or either the rows or the columns are exhausted. For those newly detected features that don’t match to any feature in the database, they are considered as potential features with unknown chemical identity.

Targeted peak detection

The targeted peak detection is conducted one feature at a time. First, a tolerance range is determined for the target feature. The m/z tolerance δ_m is determined based on the machine’s resolution, either using the nominal value or using the level found by apLCMS in its adaptive search ⁹. The retention time tolerance δ_r is determined by taking the maximum over the standard deviation of historical data times a constant of user choice, and the variation level in the current data found during peak alignment. We use x to denote m/z value, t to denote retention time, and y to denote the signal strength. Let the target m/z be m, and the target retention time be r, the target region is the rectangular region: x ∈[m − δ_m, m + δ_m], t ∈[r − δ_r, r + δ_r].

Second, the method selects a neighborhood in the profile that is much larger than the target region. Within the neighborhood, the method first iteratively split the m/z dimension into smaller sub-regions using kernel smoothing (Fig. 3). In every iteration, a smaller window size that depends on the current sub-region width is used, until a pre-specified minimum window size value is met. The iteration stops when no valley is detected in any sub-region. Then for every sub-region, the same procedure is applied in the retention time dimension. After the two iterations complete, the neighborhood is split into a number of rectangular sub-regions. Then the run filter described in ⁹ is applied to each sub-region, and peaks are called only when signals within the sub-region passes the filter. The mean m/z value weighted by intensity, Σ_i x_iy_i/Σ_iy_i, is taken as the m/z of the peak. Similarly the weighted mean retention time, Σ_i t_iy_i/Σ_iy_i, is taken as the estimate of retention time.

Third, a peak is considered to fall into the target region if its apex is in the region. Combined with the run filter, we effectively reduce false positives that result from pure noise, or the tail of a nearby peak. As illustrated in the top panel of Fig. 4, the tails of two nearby peaks fall into the target region. Because their peak apexes are outside the region, the method determines the feature is not present (Fig. 4, top panel). If one detected peak falls into the target region, then the feature is considered to be detected. If more than one peak fall into the target region, we take the apexes of the peaks, and find their Mahalanobis distance to the center of the target region, using $\left(\begin{array}{cc}{\delta }_m^2& 0\\ 0& {\delta }_r^2\end{array}\right)$ as the variance-covariance matrix. As illustrated in the lower panel of Fig. 4, both peaks #2 and #3 fall into the tolerance region. Given that peak #3 shows a smaller Mahalanobis distance to the center of the target region, it is identified as the target feature. For some features in the database, the retention time information is unavailable, because the features are derived from known metabolites that have not been observed in real data before. In such situations, the retention time tolerance range is the entire range of the experiment. If more than one peak is found within the m/z tolerance range, the one with m/z closest to the theoretical value is chosen.

Updating the database and incorporating new features

Once a new dataset is processed, the database needs to be updated to include the new information. During the previous steps, a correspondence between features in the database and those detected in the new data is already established. At the feature level, there are several scenarios.

1. The feature is in the database and not found in the new data. In this case, only the detection frequency of the feature is updated.

2. The feature is newly found and does not match to any known feature in the database. In this case, the measurements are summarized and the new feature is added to the database as a new entry.

3. The feature is in the database and detected in some of the new data. In this case, all the values need to be updated. We take retention time as an example. Assume in the database, the mean value is τ, the standard deviation is s, the feature was detected in m of M profiles previously, and we observe a new set of retention times of the feature (t1, t2,…, tn). The updated mean value is

The updated standard deviation is

Updating the range and detection frequency is straight-forward.

The m/z and log intensity values are updated in the same manner as the retention time, except when a feature is of known chemical identity, its m/z value is the theoretical value and never updated. For signal strength, there could be systematic shifts due to experiment-specific factors. This issue can be addressed by normalization using the total signal strength summed over all the detected features, before updating the signal strength information in the database.

SRM 1950 data

The Standard Reference Material (SRM) 1950 - Metabolites in Human Plasma, is made available by the National Institute of Standards and Technology (NIST). It is based on a mixture of human blood plasma. To date around 100 metabolites in the sample detectable by LC/MS have been characterized by NIST. We received the list of 94 characterized metabolites with ion forms detectable in LC/MS from Dr. Paul Rudnick’s group at NIST. For this proof-of-concept study, we constructed the known feature database using only the 94 metabolites. We analyzed the SRM 1950 sample using anion exchange chromotography combined with the Thermo Orbitrap-Velos mass spectrometer. The experiment was repeated 8 times, each in triplicate, at discrete time points spanning a month. A total of 24 spectra were generated.

We analyzed the data using three methods. The first was the untargeted feature detection using apLCMS ⁹. The second was the hybrid feature detection approach proposed in this manuscript, which is implemented in the new version of apLCMS. The third method was the untargeted feature detection using one of the leading packages XCMS ⁸. As the default setting of XCMS didn’t perform well on the data, we made an effort to tune XCMS in order to optimize its performance. Five parameters were tuned: step (values used: 0.005, 0.01, 0.05, 0.1), fwhm (values used: 15, 30, 60), mzdiff (values used: 0.05, 0.1, 0.3, 0.6), bw (values used: 15, 30, 60), and mzwid (values used: 0.005, 0.01, 0.05, 0.1). All other parameters were set at default values. The 576 possible combinations of the parameter values were all tested. We selected three parameter combinations by three different criteria: (1) The combination that yielded the maximum number of m/z matches to the 94 characterized metabolites in the sample at 10 ppm tolerance level: step=0.01, fwhm=30, mzdiff=0.05, bw=15, and mzwid=0.1. (2) The combination that yielded the maximum number of m/z matches to the 94 characterized metabolites in the sample at 5 ppm tolerance level: step=0.1, fwhm=30, mzdiff=0.6, bw=60, and mzwid=0.05. (3) The combination that yielded highest number of features detected: step=0.1, fwhm=15, mzdiff=0.01, bw=15, and mzwid=0.1. For all three methods, we required a feature to be present in at least 25% of the LC/MS profiles.

We compared the performance of the methods based on the number of features detected, the percentage of features matched to the Madison-Qingdao Metabolomics Consortium Database (MMCD) ³¹ based on mass, and the recovery of the 94 metabolites characterized by Dr. Paul Rudnick’s group at NIST (Table 1). The new hybrid approach detected the largest number of features (2115). The untargeted peak detection of apLCMS detected slightly less features (2077). The number of features detected by XCMS is much lower (maximum 1279). Secondly, we matched the features at 5 ppm level to compounds in MMCD. The allowed ion types include M, [M+H]⁺, [M+Na]⁺, [M+K]⁺, [M+NH₄]⁺. The matching was conducted using monoisotopic masses with ¹²C, ¹⁴N only. The new hybrid approach yielded the highest proportion of features matched to MMCD (39.3%). The untargeted peak detection by apLCMS followed closely (37.8%). The proportion of features matched is much lower by XCMS (maximum 30.9%). Thirdly, we considered the recovery of the 94 characterized metabolites in the sample. At 5 ppm tolerance level, the untargeted peak detection using apLCMS detected 16 of the 94 metabolites in ≥25% of the LC/MS profiles. Even lower numbers were detected by XCMS (maximum 9 out of 94). However using the new hybrid approach, 82 of the 94 features were detected in at least 25% of the profiles (Table 1). The majority of the 82 detected features were found in ≥20 of the 24 profiles.

We then manually examined the extracted ion chromatograms (EIC) of some of the features found by targeted search but not by untargeted peak detection. In fact many of such features were extremely sharp peaks (Fig. 5a). They were undetected in untargeted search because their peak characteristics (mostly small standard deviation) didn’t fit the model for general peaks. Some other features were of low ion intensity and surrounded by noise (Fig. 5b). They were undetected due to the noise filters. The results indicate that prior knowledge on known metabolites can greatly increase the detection rate of such metabolites. In this section, we only used 94 characterized metabolites for the database. In the next section, we conduct a more realistic proof-of-concept study, using a combination of HMDB ²⁹ and features from pre-existing data to build the database.

Human blood plasma data from healthy individuals

Here we present a proof-of-concept study using two previously published datasets. For the hybrid feature detection approach to work well, accumulation of feature information from large numbers of datasets of similar nature is necessary. To utilize the hybrid approach, a database needs to be built specifically for the liquid chromatography and ionization technique at hand, preferably with data generated from the same machine. In this study, we used a combination of HMDB ²⁹ and the example dataset at the apLCMS website ⁹ to initiate the database. 3286 unique m/z from HMDB and 1018 identified features from real data were incorporated into the initial database. Between the two sources, 133 of the features overlap.

We performed feature detection in the Johnson dataset which contains 30 profiles ³⁰. The dataset contains repeated measurements on 4 healthy subjects, one with 6 samples and the other 3 with 8 samples each. We required a feature to appear in at least 20% of the profiles to be considered reliable. A total of 2441 features were found in the Johnson dataset, among which 924 overlap with the Yu dataset, and 475 overlap with the HMDB database (Fig. 6a). Essentially most (90.7%) of the features found previously in the Yu dataset were also found in the Johnson dataset, showing good reproducibility of the technology, as well as the large agreement of plasma metabolites between individuals. We then categorized the newly found features into two groups: 2018 (82.3%) of the features could be directly detected without the help of the database, and 423 (17.7%) were found only through the targeted approach using the database (Fig. 6b). Among the features found only through targeted search, 274 (64.8%) were based on HMDB metabolites alone. In contrast, among the directly detectable features that match to the database, only 81 (9.5%) were matched to the [M+H]⁺ ions of HMDB metabolites.

We then studied the signal intensity and frequency of detection of the features. We compared the density plots of signal strength from features found directly and those found only by targeted search. Compared to directly detectable peaks, the peaks found only in targeted search tend to occur less in the high-strength range, and more in the low-strength range, although in the medium-strength range the two distributions are similar (Fig. 7a). Because apLCMS uses a run filter to detect features, and many medium- and low-strength peaks have missing intensities in them, it is possible that some features are missed even if their intensities are not very low. Regarding the frequency of detection among the 30 LC-MS profiles, there is a clear distinction – directly detectable features (red bars) tend to occur much more frequently than those found only by targeted search (blue bars) (Fig. 7b), which is expected. By splitting features into groups using detection frequency and generating boxplots of signal strength, we examined the data for any dependency between detection frequency and signal strength. No dependency is present in either features found directly (grey boxes) or those found only by targeted search (white boxes) (Fig. 7c).

The Johnson dataset contains repeated measures of plasma metabolites from four individuals. The directly detected features carry the information to clearly separate the individuals ³⁰. We asked the question – do the features found only by targeted search carry similar information? We conducted a principal component analysis (PCA) using only the 423 features found by targeted search. The first three PCs clearly separated the samples (Fig. 8), similar to the directly detected features ³⁰. This result indicates that targeted detection does find extra features that are biologically meaningful.

The LC/MS data is noisy. In order to identify more real features without detecting too many spurious peaks, the hybrid feature detection approach uses a simple consideration: areas corresponding to known features should receive higher confidence and thus be treated differently from other areas in terms of peak detection criteria stringency. Comparatively, if we use loose criteria all across the spectrum, we may sacrifice specificity, i.e. detect too many spurious peaks. We compared the hybrid feature detection approach with unsupervised peak detection with loosened criteria. As shown in Table 2, the unsupervised approach with loose criteria detected 4448 features, 1869 of which overlap with the 2441 features detected by the hybrid approach. Among the extra 2578 peaks the unsupervised approach with loose criteria detected, 29.4% could match at 5 ppm level to at least one compound in MMCD, while among the peaks the hybrid approach detected, 47.4% matched to MMCD. Among the 572 peaks detected by the hybrid approach and missed by the untargeted approach with loose criteria, 58.7% matched to MMCD. Although the rate of match to MMCD is not a perfect criterion, it in general reflects the quality of the features detected. These results indicate that the hybrid approach has much higher specificity in feature detection. Apparently the hybrid approach missed some true features. This was caused by the limitation of the primitive database we used in this proof-of-concept study. We only allowed [M+H]⁺ derivatives of known human metabolites in HMDB, and known features from a small dataset. With the accumulation knowledge from real data, the database will be improved and the sensitivity of the hybrid detection will increase.

In the proof of concept study, the known feature database contained 4171 entries. After processing the Johnson data, the updated database contained 5333 entries. The big increase is expected because the original database was limited – it was built with a small real dataset, as well as only [M+H]⁺ ions of HMDB metabolites. When large amounts of data are used to update the database, we expect the number of entries will increase to ten to twenty thousands. With the increase of database size, one concern is the increase of computational cost. We conducted a simulation study to examine how the method scales with the size of the database. We added artificial entries to the original database by randomly generating m/z values between 100 and 800, and increased the size of the database to 8000, 12000, 16000, 20000, and 24000 entries. We conducted the analysis at each database size setting, and recorded the computer time of the targeted feature detection (Supporting Figure 2). The LC/MS data analyzed contained 30 samples, and the computing was conducted on a machine with two 2.2 GHz quad-core Xeon CPUs and 16 Gb of RAM. The results showed that the computing time increased almost linearly up to the database size of 16000 and curved up a little with further database size increase, and remained manageable (<15 minutes for all data) at the database size of 24000 entries (Supporting Figure 2).

One difficulty in database building and metabolite matching is that some metabolites share molecular composition, hence m/z values. Without extra information, LC-MS alone in some cases cannot determine a one-to-one correspondence between found features and such metabolites. Our solution to this problem is to center the database around features, and allow each feature to have multiple possible matches to known metabolites (Supporting Figure 1). If two features share m/z and differ in retention time, they both share the same list of possible matches to known metabolites. When extra information is available, the potential matches of each feature could be reduced, and the corresponding database entries can be split. The extra information could come from LC-MS/MS or extracting retention time information from elution curves of pure chemicals using the same LC system. It is beyond the scope of the current manuscript. In the output of apLCMS, a table is provided to record the matching information between detected features in the dataset and the features in the database. Together with the feature table of each spectrum, they provide a means to trace back the measurement information (m/z, retention time & intensity) of each database entry in each sample processed. When a database entry split is necessary, such trace back can be used to update the measurement information, such as intensity and detection frequency etc.

In this proof-of-concept study, we used data from a relatively simple setup – anion exchange combined with ESI, where [M+H]⁺ ion form dominates. For other experimental setups, other ion derivatives may need to be considered in building the database, and additional steps such as compound spectra deconvolution may be necessary in the analysis of the data ²⁴.

Overall, this proof-of-concept study demonstrated that extra information is indeed recovered by using the hybrid approach of feature detection. This study used a small dataset to initiate the database. If more data is used the build the database, we can expect more features to be in the database, and hence more features found through targeted search in the new dataset. A large proportion of the features found by targeted search match to known metabolites, which suggests that using the hybrid approach can improve the number of known metabolites detected, and facilitates downstream functional analysis at the metabolic network level. To use the hybrid approach, it is critical that elution characteristics of metabolites remain the same across the datasets. This essentially requires the same column and elution gradient to be used. The most conservative practice is to keep a database for each experimental setting. Nonetheless, it is possible to use the hybrid approach across experimental settings, in which case the retention time information can no longer be used. The user can simply change the values in the retention time columns in the database to N/A, and rely on m/z alone to re-initiate the database for the new experimental setting.

Generally, the hybrid approach serves three purposes – (1) increase confidence level for lower intensity features by pooling information from historical data and known metabolites; (2) increase the number of features detected in a given dataset; (3) accumulated information about metabolite variation over large number of experiments. In many research settings, the same LC-MS procedure is repeatedly used on the same type of data. An example would be the use of LC-MS to analyze blood plasma for multiple biomarker studies. In such situations, the hybrid procedure can be repeatedly used and the database constantly updated. The accumulation of information not only helps to boost the sensitivity of feature detection, but can also help find rare features and/or features with uncommon concentration in a given set of samples. Such information can be very useful in the study of metabolic and environmental factors in diseases.