A flexible statistical model for alignment of label-free proteomics data - incorporating ion mobility and product ion information

Label-free proteomics

In a standard “bottom-up” proteomics experiment, proteins are first digested into peptides by a proteolytic enzyme. Peptides in this mixture are then physically separated by Chromatography, often Liquid Chromatography (LC). Eluting peptides are converted to gas phase ions, which are separated in a Mass Spectrometer (MS) by mass-to-charge ratio, and the relative abundance of each ion is measured by a detector. LC-MS experiments utilize a single mass analyzer, resulting in a retention time, mass-to-charge ratio, and intensity for each analyte. In LC, tandem MS experiments, or LC-MS/MS, select precursor ions are further fragmented into product ions, resulting in an additional level of information for each peptide ion. The product ions are analyzed to determine a peptide sequence, which is used to identify the parent protein. A recent variation of LC-MS/MS - Data Independent Acquisition (DIA) - generates product ions for virtually every precursor ion, providing tremendous utility for quantification and identification in a single data set. Examples of DIA include SWATH [1] and MS^E[2]. In MS^E, precursor ions enter a collision cell, rapidly alternating between high and low kinetic energy states. This “high-low switching” fragmentation enables the measurement of both precursor and product ions in a single experiment. An even more recent DIA approach to bottom-up proteomics experiments - HDMS^E - incorporates Ion Mobility (IM) spectrometry, an additional separation of peptide ions after LC, and before MS^E. IM spectrometry separates ionized peptides based on charge and three-dimensional cross-sectional area.

Label-free proteomics data processing

Several data processing steps are required to elucidate individual peptide intensities from raw label-free proteomics data. A typical data processing pipeline for a label-free proteomics experiment with multiple samples is illustrated in Figure 1. Peptide peaks must be discerned from noise, charge states determined, and isotopic distributions identified and often combined into peptide features. Further details regarding current peak detection, de-isotoping, and charge state detection methods are described in Dowsey et al. [3] and Zhang et al. [4]. The LC retention times and elution order of peptides often shift between runs. Such variations in retention time are typically called warp. The process of correcting these distortions to allow accurate matching across runs is called de-warping. Many de-warping methods exist, performing linear or non-linear (or both) corrections of two or more samples [5]. This de-warping step is either performed on raw profile data (prior to or independent of peak detection and de-isotoping), or on feature data (detected peptide features). After generation of a peptide feature set, peptide identifications are made wherever possible, and intensity measurements of both identified and unidentified peptides are grouped across runs, creating a peptide-by-sample intensity array for subsequent analyses. It should be noted that the order of data processing steps may vary within different pipelines. These data processing steps pose significant computational challenges, and are thought to be the source of much irreproducibility. This was illustrated by a recent test study by Bell et al. [6]. In the study, a sample of 20 proteins was distributed to 27 different labs, experimentally analyzed, and subjected to a variety of computationsl data-processing methods. There were significant discrepancies in reported proteins, however, all raw data was sufficient to identify all 20 proteins when centrally re-processed.

Label-free proteomics data alignment

We focus on the problems of simple data de-warping, and matching peptide intensities across multiple high-throughput proteomics runs - a combined processing step we call alignment. Our analysis emphasizes the data matching step. Accurate alignment is essential in large-scale proteomics experiments, particularly in biomarker discovery where the comparative nature of these studies require intensities of the same peptide to be compared across samples [7, 8]. In addition, accurate matching across samples can increase identifications as information can be leveraged from all individual runs [9]. The complexity of biological samples, however, poses significant computational challenges for both data alignment and peptide identification. Most samples contain tens of thousands of peptides and measurements often reflect “overlapping” peptides, or co-eluting peptides having nearly the same mass-to-charge ratio. These overlapping peptides complicate and often prevent identification. However, recent experimental advancements provide additional separation information that has not yet been leveraged in data alignment - namely comprehensive product ion information with DIA, and IM drift times with HDMS^E. Product ions have been used extensively in database matching for peptide identification [10, 11], but are not widely used in proteomics data alignment. Matching across samples is typically performed using experimentally measured and inferred characteristics of each peptide feature. Measured characteristics include precursor ion retention time, mass-to-charge ratio, and intensity. Depending on the experimental methods, some or all peptide features may have additional measured characteristics including the intensities, mass-to-charge ratios, and retention times of product ions. In DIA experiments, virtually every peptide feature has measured product ion data. In HDMS^E experiments, all precursor and product ions have a measured IM drift time as well. Inferred characteristics include charge state for nearly every peptide feature, and an amino acid sequence for some peptide features. Modern high-throughput proteomics experiments offer a great deal of information, not all of which is currently utilized in data processing - specifically in data alignment steps.

Previous alignment approaches

Matching methods utilize various aspects of the data to group peptide measurements across samples. Previously utilized characteristics include retention time, mass-to-charge ratio, intensity, and amino acid sequence. Incorporating additional data provides a higher degree of specificity when making matches [5]. The majority of existing alignment techniques utilize mass-to-charge ratio and retention time information. Such methods include SpecArray [12], AMT tag approaches [13], Xalign [14], MZmine [15], msInspect [16], XCMS [17], PETAL [18], OpenMS [19, 20], apLCMS [21], and MZmine2 [22]. Semi-supervised approaches such as PEPPeR [23] and a method by Fischer et al. [24], take advantage of existing MS/MS peptide identifications. More recent alignment methods by Tang et al. [25] and Zhang et al. [26] utilize peptide intensity along with mass-to-charge ratio and retention time. SuperHirn [27] indirectly incorporates intensity information during multiple alignments by performing pairwise alignments in a specific order based on LC-MS similarity and intensity correlations.

In addition to utilizing novel data characteristics, our method is designed to avoid common pitfalls of other approaches including static distance cutoffs, elution order assumptions, and the selection of a single reference sample. Methods should not rely on the fact that peptides always elute in the same order from the LC column as this is a false assumption and can introduce matching errors [28-30]. Multiple alignment methods often require a reference run to which all other runs are aligned. While this approach is successful for the de-warping step, choosing a single reference run for the matching step can be problematic. If the measurement variability is high in the reference sample, this can result in incorrect matches. We present a novel statistical alignment method that corrects for linear global variation, is not restricted by static distance cutoffs, and has the ability to utilize retention time, mass-to-charge ratio, peptide identifications, and previously ignored aspects of proteomics data - ion mobility drift times, and product ion data. Our method is an adapted Bayesian Dirichlet Process Gaussian Mixture Model (DPGMM) [31, 32], adding sample-specific shift and scale parameters. The proposed method is also easily extensible to incorporate additional dimensions, such as a second LC separation. We present the results of our alignment model on various datasets, comparing alignment accuracies with the inclusion of various data characteristics.

E. colilysate data

The MS^EE. coli lysate data was aligned using PEPPeR, OpenMS, and our alignment method, utilizing precursor ion mass-to-charge ratio and retention time (MZ-RT). Each alignment method was provided with the same 15 “given matches” to initialize model hyperparameters - the 15 identifications shared by all three replicates having the highest average ProteinLynx Global SERVER (PLGS) Identity^E[34] peptide score. When assessing alignment performance, we examine matches between each replicate pair - ID0821901 vs. ID0821902, ID0821901 vs. ID0821903, and ID0821902 vs. ID0821903. Correct matches occur when two identified peptides shared between the pair of technical replicates are aligned. Mismatches occur when two identified peptides with conflicting identifications are aligned. Since PEPPeR allows multiple peptides from a single sample to be present in a “cluster”, there may be more than one identification from a single replicate. In these cases, we counted a correct match if shared identifications were present in the same cluster. Similarly, we counted a mismatch if conflicting identifications were present in the same cluster. As there are varying levels of confidence in peptide identifications, we present the results of each alignment method considering identifications having a PLGS Identity^E[34] peptide score of five, six, and seven or greater. Panels A and B of Figure 2 show the recall rate and the mismatch rates, respectively, for each alignment method. As PEPPeR does not directly report match confidence, recall and mismatch rates were computed from all reported matches for each method. We see that our alignment method obtains significantly higher recall rates than OpenMS FeatureLinker and PEPPeR, for identified peptides of each confidence level. When comparing the mismatch rates, computed as the number of mismatches divided by the total matches, we see that our method obtains mismatch rates comparable with OpenMS, while PEPPeR obtains significantly higher mismatch rates, particularly at lower peptide scores. It should be noted that this mismatch rate is not a false positive rate. The total match count includes many matches which cannot be identified as correct or incorrect, as neither peptide has a putative peptide sequence. When examining the total match counts in Panel C of Figure 2, we see that our method obtains match counts comparable to PEPPeR, and significantly more matches than OpenMS.

The HDMS^EE. coli lysate data was aligned using different data characteristics to compare their utility for alignment. As with the MS^E data, each alignment was provided with the same 15 “given matches” to initialize model hyperparameters and the remaining identifications were used to assess alignment performance. We present the results of each alignment only considering identifications having a PLGS Identity^E[34] peptide score of five or greater. Results at additional peptide score thresholds are provided in Additional file 1. Comparing the HDMS^E alignments with various data combinations informs about the utility of each dimension in data alignment. We examine matches between each replicate pair - ID0822001 vs. ID0822002, ID0822001 vs. ID0822003, and ID0822002 vs. ID0822003. A match across two samples occurs when a measured peptide feature from each sample is assigned to the same latent peptide feature. Panels A and B of Figure 3 show the recall rate and the rate and number of mismatches, respectively, for each of the four alignments. We performed a series of two-sample t-tests assuming equal variance for the recall rates and mismatches to assess differences between alignments. When examining the results of the two-dimensional parent ion alignments, MZ-RT and MZ-IM, we see that the MZ-RT alignment obtains significantly higher recall rates for matches of all confidence levels. Utilizing all parent ion data characteristics obtained via HDMS^E with a three-dimensional alignment (MZ-RT-IM) results in significantly higher recall rates than the two-dimensional alignments, with a small increase in mismatches from the MZ-RT alignment, and a significant increase in mismatches from the MZ-IM alignment. When including product ion profiles, we see significant increases in recall rates, particularly with more confident matches. We see an insignificant increase in mismatches from the MZ-RT and MZ-RT-IM alignments. The alignment including product ion profiles results in much more confident matches overall. It should be noted that the results presented for the E. coli lysate analysis assume that all peptide identifications having an Identity^E[34] peptide score 5 or greater are correct. We also assume that a match of two peptides differing only by a leucine vs. isoleucine amino acid call or by amino acid order in the peptide sequences, still represents an mismatch. The resulting p-values across the range of match probability stringencies are available in Additional file 1. We also examined the match probabilities of all shared identifications - the peptides that should be matched - between the replicates, for each of the four alignments. Histograms of the match probabilities are shown in Panel C of Figure 3. We see many more confident match probabilities (near 1) of the shared identifications for the MZ-RT-IM when comparing to both two-dimensional alignments. When including product ion profiles with the MZ-RT-IM-HE alignment, we also see an increase in confident match probabilities, with a migration of all intermediate match probabilities to near-zero.

Human with E. colilysate decoy

In order to assess alignment performance without experimental identification bias, we aligned technical replicates of E. coli lysate data with a Human plasma decoy, and technical replicates of Human plasma with an E. coli lysate decoy. One technical replicate of the E. coli lysate sample was aligned with another E. coli lysate technical replicate combined in silico with a decoy Human plasma sample. To combine the samples, we append the human plasma peak list onto the peak list of one of the E. coli replicates. The Human plasma was a pooled sample from 20 individuals, run with 2 technical replicates. Similarly, one replicate of the human plasma was aligned with another plasma technical replicate combined with an E. coli lysate sample. As with the E. coli lysate data, we performed the four different alignments described in Table 1. Each alignment was provided with the same 15 “given matches” - the 15 identifications shared by the technical replicates, and having the highest average PLGS Identity^E[34] peptide score. To assess alignment performance, we determine the proportion of incorrect species alignments. This provides a false discovery rate that avoids experimental identification bias - assuming the set of shared peptides across species is negligible [35]. Figure 3, Panel D shows the correct and incorrect species match counts for each of the four alignments. We see that the MZ-RT-IM and MZ-IM alignments yield similar results, while the MZ-RT alignment obtains fewer incorrect species matches, but fewer matches overall. The MZ-RT-IM-HE alignment obtains the least incorrect species matches, and many more correct species matches at reasonably confident match levels. The results of the inverse decoy analysis (aligning technical replicates of E. coli lysate with a Human plasma decoy) are available in Additional file 1: Figure S7.

Hepatitis-C and osteoarthritis data

To illustrate the utility of aligning datasets obtained from two different tissues, we aligned MS^E serum samples of Hepatitis-C patients with MS^E human urine samples from an Osteoarthritis cohort using peptide ion mass-to-charge ratio and retention time. Peptide identifications from the urine samples were carried over to the serum via alignment. We examined the 500 peptides having the most significant differential expression with a phenotype of interest - Hepatitis-C treatment response. Significance was assessed with two-sample t-tests assuming equal variance, on peptide intensities that were log-transformed, mean-centered by sample, and standardized by peptide. We explored the inferred identifications that were carried over via alignment from urine, in order to identify potential biomarkers. Specifically, we looked at the previously unidentified peptides exhibiting significant differential expression having an inferred identification from the alignment results. We analyzed these inferred identifications using GATHER [36] and DAVID [37] to search for functional and pathway enrichment. The lists of proteins from these inferred identifications and their GATHER and DAVID analysis results are available in the Additional file 1. The inferred identifications of peptides exhibiting differential expression for Hepatitis-C treatment response totaled 42 corresponding proteins. These proteins were significantly enriched for defense response (GO0006952 - 15 proteins), immune response (GO0006955 - 15 proteins), and response to biotic stimulus (GO0009607 - 14 proteins), having GATHER p-values 1.56e-9, 3.76e-9, and 9.81e-9, respectively. These functional annotations are very much in accordance with what one would expect with a response to viral infection, suggesting that useful identifications were obtained with the alignment. In addition, the genes encoding 6 of these proteins were located at 14q32, indicating significant chromosome location enrichment with GATHER p-value 1.27e-5. This is the location of the immunoglobulin heavy locus - a region containing genes encoding the heavy chains of antibodies. Differential expression of genes in this region is also in accordance with what one would expect for HCV treatment response.

Alignment model

We present a statistical model for the alignment of label-free proteomics data to match peptide features across multiple samples after peak-detection and de-isotoping. Unlike any existing proteomics alignment method, our model has the ability to utilize ion mobility drift time from HDMS^E experiments, and product ion spectra from traditional LC-MS/MS Data Dependent Aquisition (DDA), or DIA (MS^E or HDMS^E) experiments, along with the typical parent ion mass-to-charge ratio and retention time - increasing the individuality of each peptide feature and providing a better alignment. At the time of publication, no open-source proteomics file format was capable of storing ion mobility separation data. In order to allow incorporation of this data into our alignment method, we wrote a small data-processing script to read Waters ‘spectrum.xml’ and ‘finalfrag.csv’ files into a Matlab data frame. The Matlab data frame format, a sample configuation file, a sample queue submission file, and the data processing script are available in Additional files 1, 3, 4 and 5, respectively, and can be easily adapted to incorporate any additional separation dimensions similar to ion mobility drift times and liquid chromatography retention times, including retention times from multidimensional LC.

Peptide-level model

A sample-specific linear shift (η _d ) and scale (β _d ) is used for de-warping. In the formulas that follow, samples or datasets are indexed by d, individual measured peptides within a sample are indexed by i, and latent (theoretical) peptides are indexed by j. As shown in Equation 1, we assume that a measured peptide feature, x_d,i, is a shifted and scaled noisy measurement of a “true” peptide feature, $z_{c_{d,i}}$. Let c_d,i be an indicator variable for the latent peptide assignment of measurement x_d,i, taking on values j=1…J where J is unbounded.

$\begin{array}{l}{x}_{d,i}={\eta }_d+z_{c_{d,i}}{\beta }_d+{\epsilon }_{d,i}\end{array}$

(1)

$\begin{array}{l}{\epsilon }_{d,i}\sim \mathit{\text{Normal}}(0,\Sigma )\end{array}$

(2)

$\begin{array}{l}{z}_{c_{d,i}}\sim \mathit{\text{Normal}}({\mu }_{c_{d,i}},\sigma )\end{array}$

(3)

Where ε_d,i are the residuals between the measured values (x_d,i) and the shifted and scaled latent values (${\eta }_d+z_{c_{d,i}}{\beta }_d$), having a multivariate normal distribution - Equation 2. Let Σ be the covariance of these residuals, and let each latent peptide $z_{c_{d,i}}$ be a draw from DPGMM mixture component c_d,i=1…J having mean ${\mu }_{c_{d,i}}$ and covariance σ, as shown in Equation 3. Note that we make the simplifying assumption of shared covariance across all latent peptides. This would suggest that the measured values of each latent peptide show the same variation across the entire mass-to-charge ratio, retention time, and drift time range. We acknowledge that this is not likely the case, although we find this assumption works well in practice. Conjugate priors are used for all model parameters as follows:

$\begin{array}{l}{\eta }_d\sim \mathit{\text{Normal}}(a_d,b_d)\end{array}$

(4)

$\begin{array}{l}{\beta }_d\sim \mathit{\text{Normal}}(e_d,f_d)\end{array}$

(5)

$\begin{array}{l}\Sigma \sim \mathit{\text{Inverse}}-\mathit{\text{Wishart}}(s,t)\end{array}$

(6)

$\begin{array}{l}{\mu }_j\sim \mathit{\text{Normal}}(\lambda ,r)\end{array}$

(7)

$\begin{array}{l}\sigma \sim \mathit{\text{Inverse}}-\mathit{\text{Wishart}}(g,h)\end{array}$

(8)

Normal priors are assigned to the shift and scale parameters as shown in Equations 4 and 5. The seed matches, for each peptide feature are averaged to generate a list of “implied-identified peptides”. Robust fit linear regression is performed for each dataset using the “implied-identified peptides” as predictors, and the measured identified peptides as response. The resulting intercept is taken as the mean hyperparameter in the shift prior distribution (a _d ). Similarly, the coefficient is taken as the mean hyperparameter in the scale prior distribution (e _d ). The variance parameters on the shift and scale priors (b _d and f _d ) are set tightly to the variance of the regression estimate. This allows an optimal solution to be reached as latent peptides are updated and added, while reducing shift and scale identifiability issues. Both the match covariance (Σ), and latent peptide covariance (σ) matrices are given conjugate inverse-Wishart priors as shown in Equations 6 and 8. The residuals of the shifted and scaled identified peptide measurements, and their respective “implied-identified peptides” are used set hyperparameters. The degrees of freedom parameters (h and t) are set to the number of identified matches minus one, and the inverse-scale matrix is set to the sum of squared residuals. The mean of each latent peptide (μ _j ) is given a conjugate normal prior with as shown in Equation 7. The prior mean (λ) is set to the empirical mean of all measured peptide features in all datasets, and the prior covariance (r) is set to the sum of squared differences between this empirical mean and all measured data. We express the likelihood of x_d,i as follows:

$\begin{array}{ll}P(x_{d,i}\mid & \eta ,\beta ,\Sigma ,z_1\mathrm{...}{z}_J,c_{d,i})=\mathit{\text{Normal}}(x_{d,i}\mid {\eta }_d+z_{c_{d,i}}{\beta }_d,\Sigma )\end{array}$

(9)

Where z₁...z _J indicates all existing latent peptides. We may also integrate out $z_{c_{d,i}}$ and re-express the likelihood as:

$\begin{array}{ll}P(x_{d,i}\mid & \eta ,\beta ,\Sigma ,{\mu }_1\mathrm{...}{\mu }_J,\sigma ,c_{d,i})=\mathit{\text{Normal}}(x_{d,i}\mid A_{c_{d,i}},B_d)\end{array}$

(10)

$\begin{array}{ll}{A}_{c_{d,i}}& ={\eta }_d+{\mu }_{c_{d,i}}{\beta }_d\\ B_d& ={\beta }_d^T\sigma {\beta }_d+\Sigma \end{array}$

The prior probability of an observation, x_d,i, being assigned to latent peptide component j given all other assignment indicators, c_−(d,i), is given in Equation 11. The notation c_−(d,i) refers to all component indicators from all features in all datasets, except d,i. Similarly the prior probability of observation x_d,i being assigned to a new latent peptide component is shown in Equation 12.

$\begin{array}{l}p(c_{d,i}=j\mid c_{-(d,i)},\alpha )=\frac{n_{-(d,i),j}}{N-1+\alpha }\times I(!\exists i,c_{d,-i}=c_{d,i})\end{array}$

(11)

$\begin{array}{l}p(c_{d,i}\ne c_{-(d,i)}\mid c_{-(d,i)},\alpha )=\frac{\alpha }{N-1+\alpha }\end{array}$

(12)

Let α be the DPGMM concentration parameter, N the total number of observed peptide features across all samples, and n_−(d,i) the number of observed peptide features other than x_d,i assigned to latent peptide j. A constraint is imposed such that only one measurement per dataset may be assigned to a given latent peptide feature. The concentration parameter of the Dirichlet Process (α) is set to the number of peptide feature observations across all samples being aligned. Latent peptide feature assignments are updated from their full conditional posterior distributions, as shown in Equation 13.

$\begin{array}{ll}P(c_{d,i}=j\mid -)& \propto \frac{\alpha }{N-1+\alpha }\times \int P(x_{d,i}\mid -)\\ \times P({\mu }_j\mid \lambda ,r)d{\mu }_j\\ +\sum _j\frac{n_{-(d,i),j}}{N-1+\alpha }\mathit{\text{Normal}}(x_{d,i}\mid A_{c_{d,i}},B_d)\\ \times I(!\exists i,c_{d,-i}=c_{d,i})\end{array}$

(13)

The integral above is tractable due to the shared covariance across latent peptide components. Further details and full conditional distributions are available in Additional file 1. Panel A of Figure 4 shows a plate diagram of the peptide level alignment model.

Production model extension

To incorporate production data, we select up to the 50 most intense productions for each peptide feature measurement, x_d,i. We then generate a K-dimensional production intensity profile for each x_d,i. Each position, $y_{d,i_k}$, in the product ion intensity profile, y_d,i, is computed as:

$\begin{array}{l}{y}_{d,i_k}=\frac{\sum _p{\Omega }_p\times I(M_p\le B_k)}{\sum _p{\Omega }_p}\end{array}$

(14)

where k=1…K, p=1…50, Ω is a 50-dimensional vector of intensities, M is a 50-dimensional vector of product ion mass-to-charge ratios, and B is a K-dimensional vector of product ion profile mass-to-charge ratio bin upper limits. All values in y_d,i sum to one. The mass-to-charge ratio ranges, or bins (B _k ), are determined at the initialization of the alignment, such that bin boundaries fall on mass-to-charge ratio “deserts”. See the Parellelization section and Figure 5 for a detailed description of bin boundary determination. In the experiments described here, we set K = 250 to ensure most product ions would be assigned to their own bin in the product ion profile, and to avoid additional computational complexity. Each existing latent peptide feature is given a K-dimensional product ion profile (w _j for latent peptide feature z _j ). To assess the similarity of a measured product ion profile, y_d,i, and a latent production profile, $w_{c_{d,i}}$, we introduce a similarity score, ψ, which is computed as the sum of squared differences of the two product ion intensity profiles, and is assumed to have an exponential distribution to encourage distances close to zero, as shown in Equations 15 and 16.

$\begin{array}{l}{\psi }_{d,i}={(y_{\mathit{\text{di}}}-w_{c_{d,i}})}^T(y_{d,i}-w_{c_{d,i}})\end{array}$

(15)

$\begin{array}{l}{\psi }_{d,i}\sim \mathit{\text{Exponential}}\left(\gamma \right)\end{array}$

(16)

We assign a conjugate gamma prior to the rate parameter, as shown in Equation 17. The hyperparameters for profile scores are set to one.

$\begin{array}{l}\gamma \sim \mathit{\text{Gamma}}(a_0,b_0)\end{array}$

(17)

At each iteration of the MCMC, the product ion profile, w _j of an existing latent peptide is updated empirically - the product ion profile is set to the average of the measured product profiles assigned to that latent peptide. The latent product ion profile, w₀, of a new latent peptide (one that currently does not exist) is a blank profile - a uniform vector of size K with each element having value 1/K. Combining the product ion model with the peptide-level model, we have the following likelihood (Equation 18) and conditional posterior (Equation 19):

$\begin{array}{ll}P(x_{d,i},y_{d,i}\mid -)& =\mathit{\text{Normal}}(x_{d,i}\mid A_{c_{d,i}},B_d)\\ \times \mathit{\text{Exponential}}({\psi }_{d,i}\mid \gamma )\end{array}$

(18)

$\begin{array}{ll}P(c_{d,i}=j\mid -)\propto \frac{\alpha }{N-1+\alpha }& \times \int P(x_{d,i}\mid -)\\ \times P({\mu }_j\mid \lambda ,r)d{\mu }_j\\ \times \mathit{\text{Exponential}}({(y_{\mathit{\text{di}}}-w_0)}^T\\ \times (y_{d,i}-w_0)\mid \gamma )\\ +\sum _j\frac{n_{-(d,i),j}}{N-1+\alpha }\\ \times \mathit{\text{Normal}}(x_{d,i}\mid A_{c_{d,i}},B_d)\\ \times I(!\exists i,c_{d,-i}=c_{d,i})\\ \times \mathit{\text{Exponential}}({\psi }_{d,i}\mid \gamma )\end{array}$

(19)

Further details and full conditional distributions are available in Additional file 1. We explored additional values of K, as well as implementations of different product ion models, the results and discussions of which can also be found in Additional file 1.

Data

All data used in this analysis was obtained under MS^E and/or HDMS^E conditions (SYNAPT HDMS G2, Waters), and subject to Waters ProteinLynx Global SERVER (PLGS) processing. We utilize peptide features that have already been subject to peak detection, de-isotoping, charge state determination, and tentative identification (although not all identifications are utilized for alignment). All samples were separated by 1D nanoscale capillary ultraperformance Liquid Chromatography in a 90-minute gradient using a 5-40% acetonitrile/water (0.1% formic acid in each).

Three HCV cohorts were utilized in the alignment of serum samples from HCV patients to urine samples from OA patients. The first cohort included 47 patients ages 5 to 18 years from a clinical trial for HCV treatment [41]. The two additional HCV cohorts (n = 41,55) were selected from the Duke Hepatology Clinical Research (DHCR) database [42]. The pediatric clinical trial study was approved by the institutional review boards of the participating sites. Written informed consent was provided by all parents or guardians, and written assent was provided by all participants over 12 years of age. All patients present in the DHCR database cohorts, as well as all OA patients, provided written informed consent, and all study procedures were approved by the Duke University Institutional Review Board.

Analysis

All alignments were performed using Matlab on the Duke Shared Computing Resource, a cluster of Intel x86 compute notes running Linux. Each alignment was partitioned into a maximum of 250 splits, and each alignment partition was run on a single node with at least 8GB of memory. Our method does require considerable computation time - approximately 2 to 6 hours for the E. coli Lysate and HCV-OA alignments not utilizing product ions, and approximately 20-24 hours for the E. coli Lysate alignments utilizing product ions. Times vary by the number and size of datasets being aligned.