Bayesian Mass Spectra Peak Alignment from Mass Charge Ratios

Abstract

Proteomics studies based on mass spectrometry (MS) are gaining popular applications in biomedical research for protein identification/quantification and biomarker discovery, especially for potential early diagnosis and prognosis of severe disease before the occurrence of symptoms. However, MS data collected using current technologies are very noisy and appropriate data preprocessing is critical for successful applications of MS-based approaches. Among various data preprocessing steps, peak alignment from multiple spectra based on detected peak sample locations presents special statistical challenges when effective experimental calibration is not feasible due to relatively large peak location variation. To avoid intensive tuning parameter optimization, we propose a simple novel Bayesian algorithm “random grafting-pruning Markov chain Monte Carlo (RGPMCMC)” that can be applied to global MS peak alignment and to follow certain modelbased sample classification criterion for using aligned peaks to classify spectrum samples. The usefulness of our approach is demonstrated through simulation study by making extensive comparison with other algorithms in the literature. Its application to an ovarian cancer MALDI-MS data set achieves a smaller 10-fold cross validation error rate than other current large scale methodologies.

1. Introduction

Genomics and proteomics technologies offer much promise in our understanding of fundamental biological processes by allowing us to simultaneously monitor the expression levels of tens of thousands of genes and proteins. Since proteins are the basic functioning units in the cells, there is a great interest to characterize individual molecular profiles based on proteomics for reliable biomarker discovery and effective disease diagnosis/prognosis/treatment. In proteomics research, mass spectrometry (MS) is the most widely used instrument to allow for the mass measurement of molecules, where a mass spectrometer determines chemical compounds’ molecular weight by ionizing, separating, and measuring molecular ions according to their mass-to-charge ratio (m/z: unit Da) and a mass spectrum is the standard data output for analysis and interpretation (Link et al. 1999), where the x-axis represents m/z value (Da) and the y-axis represents intensity (enrichment of particles with certain m/z). Recently, Yu et al. (2005) reviewed current approaches on the extraction of the most relevant information from the raw mass spectra to identify disease biomarkers. A standard MS data analysis usually involves background noise removal, smoothing, intensity normalization, peak identification and alignment, and biomarker identification. The false positive and false negative peaks may exist in taking local maxima as peaks (Coombes et al. 2003) and peak location variation may be due to differences in sample preparation, chemical noise, co-crystallization, deposition of the matrix-sample onto the target, and laser position on the target among others. Several statistical methods have been proposed to reduce the background noise (Coombes et al. 2003; Satten et al. 2004; Coombes et al. 2005a and Randolph and Yasui, 2006). Morris et al. (2005) applied translation-invariant wavelet transformations to the raw spectra and performed peak detection using the mean spectrum derived from a group of spectra, where they assumed calibration can be done in advance experimentally or by interpolation to make common peaks stay closely together, and only false negatives are possible. However, Coombes et al. (2005b) pointed out that, peak location variation, caused by a spread of initial particle velocities at the starting end of mass spectrometer tube, makes calibration more difficult. Compared to background noise removal, peak identification and alignment is more challenging and critical by providing the links of underlying peptides across all spectra. However, existing algorithms for peak alignment are mostly ad hoc and based on heuristic arguments (Nielsen et al. 1998; Johnson et al. 2003; Torgrip et al. 2003; Eilers, 2004; Tibshirani et al. 2004 and Randolph and Yasui 2006), where some parameters need to be optimized empirically and/or subjectively. In this paper, we work on the detected peak samples (possibly with associated intensities) coming from certain protocol. In view of hundreds of features (peaks) in each spectrum, due to complex chemical and physical mechanisms undergoing the mass spectrometry, throughout this paper, we assume that substantial individual peak location variation is existent, say up to one half of the interval between neighboring peaks, thus calibration may be in lack of power. Moreover, we consider statistically false positives and negatives to make our model more accountable. Overall, we will assume that each set of potential peak samples corresponding to the true peak follows an individual composite distribution regulated by true peak location, peak sample location variation, false positive and false negative rates. An effective simple Bayesian MCMC algorithm is proposed to do peak alignment (biomarker identification) and downstream sample classification, where all individual sets of parameters for each true potential peak are able to be estimated in a universal modeling framework without intensive tuning parameter optimization. Our algorithm is much computationally simpler than other available MCMC algorithms which could be applied to MS alignment while retaining competing performance.

The rest of this article is organized as follows: Section 2 introduces Bayesian dimension matching problem and proposes our simplified approach; Section 3 summarizes simulation results to demonstrate the efficiency and reliability of our algorithm; Section 4 illustrates the applications of our approach to a real MS data set; Section 5 considers joint analysis of peak alignment and sample classification; Section 6 concludes with discussions; and some technical details are given in the appendix.

2. Methods

2.1. Dimension-Matching Statistical Model

We now develop a simple novel MCMC algorithm, random grafting-pruning Markov chain Monte Carlo (RGPMCMC) which can be applied to MS peak alignment from mass charge ratio information. Within Bayesian framework, we are often interested in the posterior distribution of the dimension-varying parameter θ. The prior distribution π(θ)is represented as ${\sum }_{k\in N}\text{π}\left({\text{θ}}_K|K\right)\text{π(K)}$ , where N is positive integer set and $\text{π}\left({\text{θ}}_k|K\right)$ is the individual prior distribution within the K-dimensional space and π(K)is the mixture probability for dimension K. Since $\text{π(}K,{\text{θ}}_k\text{)=}\hspace{0.17em}\text{π(}K\text{)π(θ|K)}$ , the posterior distribution of (K,θ _K )

$$\begin{array}{l}\pi \left(K,{\theta }_k|X\right)=\\ \frac{f\left(X|{\theta }_k\right)\pi ({\theta }_k)(K)}{{\sum }_{K\in K}({\int }_{{\Theta }_k}f(X|{\theta }_k)\pi ({\theta }_k)d{\theta }_k)\pi (K)},\end{array}$$

(1)

where the denominator is the normalization constant not needed for posterior sampling. Change point model usually involves dimension-matching in two typical cases: one single ordered series where change points are taken as those separating successive discrete points, and multiple ordered series where change points correspond to physical locations in the continuous space. For the discrete case, the partitioning and wrapping up the segments leads to an exponentially increasing computational cost as the model space grows (Denison et al. 2002), and frequentist's approaches only either work on very few change points or special algorithms (Guan, 2004; Olshen and Venkatraman, 2004). Recently, Fearnhead (2005) proposed an exact non-MCMC sampler by recursive partition, and Loschi, Cruz and Arellano-Valle (2005) introduced the product partition model based Bayesian algorithm which stems from Yao (1984) and Barry and Hartigan (1992, 1993). For the continuous case, the readers are referred to sampling-based algorithms: the reversible jump MCMC (RJMCMC) by Green (1995), the birth-and-death process MCMC (BDMCMC) and continuous time process MCMC (CTMCMC) by Stephens (2000a, 2000b), Bayesian cluster detection in maps (Knorr-Held and Raßer, 2000) and others. Cappέ, Robert and Rydέn (2003) showed that, the acceptance probability of the usual MCMC methods is replaced by differential holding times in BDMCMC, and RJMCMC converges to a limiting continuous time birth-and-death process on an appropriate rescaling of time. They also demonstrated that, RJMCMC and CTMCMC have similar computational performance while the latter demands expensive death rate computation. Obviously, discretization is neither always suitable nor efficient for sophisticated change point identification in the continuous space. The present work introduces a simple MCMC algorithm in the context of multiple MS peak alignment by considering all uncertainties including peak number, peak locations, peak sample location variations, false negative and false positive rates. Neither the error-prone Jacobian terms in RJMCMC, the intensive death rate calculation in CTMCMC, nor the computationally expensive recursive partition in other algorithms is needed by our method.

Before starting with our statistical model, we assume local maxima (discrete locations) have been detected as peak samples for each raw mass spectrum. In MS peak alignment, peak sample location variation may linearly depend on m/z magnitude (Yasui et al. 2003), and log-transformation of m/z achieves peak sample location variation homogeneity, this observation will justify the identical prior specification for peak sample location variations across true peaks (details later). For notational simplicity, we use m/z instead of log (m/z) in the following discussions and assume that the m/z domain is [(m/z)_min, (m/z)_max]. The underlying K-dimensional true peak location vector is ${\tilde{S}}_k$ = (s₁, s₂, …, s_K), where the peak locations are generally separated by at least a distance threshold, say no less than d. The data to be analyzed from multiple spectra are the detected peak sample locations y_ij (j = 1, …, n_i, i = 1, …, I), where i is the index for spectra, j is the peak sample index within each spectrum (with increasing m/z), and n_i is the number of detected peak samples in the i-th spectrum. The peak samples are assumed to be normally distributed around their true peak s_k with standard deviation σ_k (k = 1, …, K) for locations. For the putative true peak set ${\tilde{S}}_k$ = (s₁, s₂, …, s_K), each detected peak sample j in i-th spectrum is assigned to its nearest putative true peak among ${\tilde{S}}_k$ , say n_y (i, j) (with possible multiple assignments to the same putative peak from a given spectrum). For certain true peak k, (1) when there is no peak sample assignment to it from a given spectrum, we consider it as a false negative case for this spectrum with probability fn_k; (2) when there are multiple peak sample assignments to it from a given spectrum, we consider it as a false positive case for this spectrum with probability fp_k; (3) otherwise we consider it as a non-false positive or negative case for this spectrum with probability 1 – fn_k – fp_k. Our prior hypothesis is that, each true peak shows these three types of cases proportionally under the homogeneity assumption for the spectra and identical peak sample detection protocol for each local m/z region, say for each true peak, on average 90% spectra will contribute single peak sample, 5% will contribute multiple peak samples and 5% will contribute no peak sample to the putative true peak. For these three cases, we assume an independent trinomial distribution Tri(I; fn_k, fp_k, 1 – fn_k – fp_k) for each true peak k in the context of numbers of grabbed peak samples from I spectra. Now we restate the following notations for model set-up: Y stands for the peak sample locations for all I spectra, which need not be a I-row matrix since the numbers of peak samples are not necessarily equal because of false negatives and/or false positives; ${\tilde{S}}_k$ is a K-dimensional vector of putative true peak locations (m/z's); ${\tilde{\sigma }}_k$ is a K-dimensional vector of peak sample location variations at putative true peaks; $\tilde{f}\hspace{0.17em}{n}_{\hspace{0.17em}k}$ and $\tilde{f}\hspace{0.17em}{p}_{\hspace{0.17em}k}$ are K-dimensional vectors of false negative and false positive rates at putative true peaks; n_y (i, j) is putative true peak assignment to peak sample j in the i-th spectrum; n_fn,k is the number of spectra without peak samples assigned to putative true peak k; n_fp,k is the number of spectra with multiple peak samples assigned to putative true peak k, and $n_{\overline{fnp,k}}$ is the number of spectra with single peak sample assigned to putative true peak k (Figure 1). The likelihood for the observed peak samples across all spectra is

Figure 1

Mass Spectrum Biomarker Model. (The left column is spectrum index, the vertical lines are putative true peaks, the horizontal circle lines are peak samples of each spectrum. The lower left dash rectangle represents a false negative case and the upper right dash rectangle represents a false positive case.)

$$\begin{array}{c}f\left(Y|K,{\tilde{S}}_K,{\tilde{\sigma }}_K,\tilde{f}{n}_K,\tilde{f}{n}_K,\tilde{f}{p}_K\right)\\ =\underset{i=1}{\overset{I}{\Pi }}\underset{1\le j\le n_j}{\Pi }\varphi (y_{ij}|s_{n_y(i,j)},{\sigma }_{n_y(i,j)})\\ \times \underset{k=1}{\overset{K}{\Pi }}f{n}_k^{n_{fn,k}}f{p}_k^{n_{fp,k}}{\left(1-f{n}_k-f{p}_k\right)}^{n_{\overline{nfnp,k}}}.\end{array}$$

(2)

We now describe the motivation for such a likelihood function which is crucial for Bayesian inference. The matching by minimum distance $(s_{n_y(i,j)})$ is to capture important clustering information. We believe that, this objective construction is reasonable compared to some alternatives, e.g. non-model based clustering algorithms, or piecewise binning method where each bin is assumed to hold exactly those peak samples for the individual putative true peak. The other part incorporating false positives and false negatives is mainly a technical statistical consideration, since it is very hard to accurately claim which peak sample is a true false positive or false negative. To sum up, in spite of substantial measurement errors underlying high-throughput mass spectra, we pursue a reasonable eclectic theme to tackle biomarker profile estimation. For notational convenience, “peak” represents putative true peak for the biomarker profile hereafter, other than detected peak sample from each spectrum. On the prior part, we assume that K follows a truncated Poisson or discrete uniform distribution on {K_min, …, K_max}. As Green (1995) suggested, the peak locations are taken as even-numbered order statistics from 2K + 1 points uniformly distributed on an L-length interval [(m/z)_min, (m/z)_max] (for convenience, s₀ = (m/z)_min, s_K+1 = (m/z)_max) to avoid too many short steps, which has density ${\Pi }_{k=1}^{K+1}(s_k-s_{k-1})/L^{2K+1}$ as suggested by Green (1995). To make use of conjugate prior, π(σ ² ) is taken as Inverse-Gamma (ν,η) density $({\sigma }^2){-}^{(v+1)}{e}^{-1/({\sigma }^2\eta )}{\eta }^{-v}/\Gamma (v);$ the joint prior distribution for (fn, fp, – fn – fp) is a 3-dimensional Dirichlet distribution with density $D\left(fn,fp,1-fn-fp|{\alpha }_1,{\alpha }_2,{\alpha }_3\right)={(fn)}^{{\alpha }_1-1}{(fp)}^{{\alpha }_2-1}(1-fn-2)fp{)}^{{\alpha }_3-1}\Gamma ({\alpha }_1)\Gamma ({\alpha }_2)\Gamma ({\alpha }_3)/\Gamma ({\alpha }_1+{\alpha }_2+{\alpha }_3)$ . The posterior distribution is

$$\begin{array}{cc}f(K,{\overline{S}}_k,{\tilde{\sigma }}_K,f{n}_K,f{p}_K|Y){\Pi }_{i=1}^I{\Pi }_{1\le j\le n_j}\varphi (y_{ij};s_{n_y(i,j)},{\sigma }_{n_y(i,j)})& (\text{Normal}\hspace{0.17em}\text{Dist}.)\\ \times \hspace{0.17em}{\Pi }_{k=1}^K\hspace{0.17em}f{n}_k^{n_{fn,k}}{(1-f{n}_k-f{p}_k)}^{n_{\overline{fnp,k}}}& (\text{trinomial}\hspace{0.17em}\text{Dist}.)\\ \times \hspace{0.17em}\pi (K,S_k){\Pi }_{k=1}^K\hspace{0.17em}[D(f{n}_k,f{p}_k,1-f{n}_k-f{p}_k|{\alpha }_1,{\alpha }_2,{\alpha }_3)\hspace{0.17em}\times \hspace{0.17em}{\Pi }_{k=1}^{K}IG({\sigma }_k^2|v,\eta )]& (\text{prior}\hspace{0.17em}\text{Dist}.)\end{array}$$

(3)

2.2. Random Grafting-Pruning Markov Chain Monte Carlo(Rgpmcmc)

Motivated by Green (1995), our algorithm is based on a redesigned universal “naively informative” parameter proposal involving peaks and the other parameters concurrently without the need for Jacobian terms. We also propose four move types (+, –, H, S), where “+” means peak birth proposal, “–” means peak death proposal, “H” means parameter (σ ² , fn and fp) proposal excluding peaks and “S” means peak location mutation with peak number unchanged. We specify (+, –, H, S) probabilities as $(\text{π}(+),\text{π}(-),\text{π}(H),\text{π}(S))$ .

Parameter Sampling Process

1. 1)

   First we choose one of these four move types based on move type probabilities (π(+), π(–), π(H), π(S)), where π(+) = π(–).

   2)

   For the “+” move type, we describe the parameter proposal process: If K_old = K_max, we go to 1) since the upper threshold is reached; if K_old < K_max, we randomly sample one of the K_old + 1 intervals formed by current K_old peaks, say (s_j, s_j+1), with equal probability 1/(K_old + 1). We may assign j + 1 to this new peak index ^* and the following indexes increase by one accordingly. Within this sampled interval, we propose ( $(s_{*},{\sigma }_{*}^2,f{n}_{*}$ and fp_*) for peak candidate ^* as follows:

   • i.

     True peak location proposal for peak birth:

     $$s_{*}=g_{s_{*}}(U_1;s_j,s_{j+1})\in \hspace{0.17em}(s_j,s_{j+1}),$$

     where g_s*(U₁; s_j, s_j+1) is a one-to-one mapping from random variable U₁ to peak location s_* given s_j and s_j+1. We take $g_{{\sigma }_{*}^2}(U_2;{\sigma }_j^2,{\sigma }_{j+1}^2)$ to be $(s_j+s_{j+1}{g}_1(U_1))/(1+g_1(U_1^{*})),$ or $(s_{*}-s_j)/(s_{j+1}-s_{*})=g_1(U_1)$ , where g₁(·) is any monotonic function with domain [0, 1] and range [0, ∞), and U₁ ~ U [0, 1]. It can be seen that, s_* is a monotonically increasing function of U₁. We simply use g₁(u) = u/(1 – u), thus s_* = s_j + (s_j+1 – s_j) U₁, a uniform random variable ∊ (s_j, s_j+1).

   • ii.

     Peak sample location variance proposal for peak birth:

     $${\sigma }_{*}^2=g_{{\sigma }_{*}^2}(U_2;{\sigma }_j^2,{\sigma }_{j+1}^2),$$

     where $g_{{\sigma }_{*}^2}(U_2;{\sigma }_j^2,{\sigma }_{j+1}^2)$ is a one-to-one mapping from random variable U₂ to peak sample location variance $g_{*}^2$ given ${\sigma }_j^2$ and ${\sigma }_{j+1}^2$ . We take $g_{{\sigma }_{*}^2}(U_2;{\sigma }_j^2,{\sigma }_{j+1}^2)$ to be ${({\sigma }_j^2{\sigma }_{j+1}^2)}^{1/2}{g}_2(U_2)$ or $({\sigma }_{*}/{\sigma }_j)/({\sigma }_{j+1}/{\sigma }_{*})=g_2(U_2)$ , where g₂(·) is any monotonic function with domain [0, 1] and range [0, ∞), and U₂ ~ U [0, 1]. It can be seen that, ${\sigma }_2^{*}$ is a monotonically increasing function of U₂. We simply use g₂(u) = u/(1 – u), thus ${\sigma }_{*}^2={({\sigma }_j^2{\sigma }_{j+1}^2)}^{1/2}{U}_2/(1-U_2)$ .

   • iii.

     Peak sample false negative and false positive rate proposal for peak birth:

     $$(f{n}_{*},f{p}_{*})={\tilde{g}}_{f{n}_{*},f{p}_{*}}(U_3,U_4;f{n}_j,f{p}_j,f{n}_{j+1},f{p}_{j+1}),$$

     where ${\tilde{g}}_{f{n}_{*},f{p}_{*}}(U_3,U_4;f{n}_j,f{p}_j,f{n}_{j+1},f{p}_{j+1})$ , is a one-to-one mapping from (U₃, U₄) to (fn_*,fp_*) given (fn_j, fp_j) and (fn_j+1, fp_j+1). Specifically, for peak ^*, we use O^* to represent “false negative or positive” odds and R_* to represent “false negative vs. positive” ratio, i.e.

     $$O_{*}=(f{n}_{*}+f{p}_{*})/(1-f{n}_{*}-f{p}_{*})\hspace{0.17em}\text{and}\hspace{0.17em}{R}_{*}=f{n}_{*}/f{p}_{*},$$

     (4)

     the false negative and false positive rate proposal is realized in two sequential steps:

   • 1)

     O_* proposal:

     $$O_{*}={(O_j{O}_{j+1})}^{1/2}{g}_3(U_3),$$

     (5)

     i.e. (O_*/O)^1/2/(O_j+1/O_*)^1/2 = g₃(U₃), where g₃(·) is any monotonic function with domain [0, 1] and range [0, ∞), and U₃ ~U[0, 1]. It can be seen that, O^* is a monotonically increasing function of U₃, we simply use g₃(u) = u/(1 – u).

   • 2)

     R_* proposal:

     $$R_{*}={(R_j{R}_{j+1})}^{1/2}{g}_4(U_4),$$

     (6)

     i.e. (R_*/R_j)^1/2/(R_j+1/R_*)^1/2 = g₄(U₄), where g₄(·) is any monotonic function with domain [0, 1] and range [0, ∞), and U₄ ~ U[0, 1]. It can be seen that, R_* is a monotonically increasing function of U₄, we simply use g₄(u) = u/(1 –u).

     Note that the constraint 0 ≤ fn_* + fp_* ≤ 1 holds under this proposal. The ${\tilde{g}}_{f{n}_{*},f{p}_{*}}(U_3,U_4;f{n}_j,f{p}_j,f{n}_{j+1},f{p}_{j+1})$ function is a combination of O_* proposal and R_* proposal in this case. fn_* and fp_* are jointly proposed to meet the constraint. The Jacobian of transforming (fn_j,fp_j,fn_j+1,fp_j+1,u₃, u₄) into ((fn_j,fp_j,fn_*,fp^*,fn_j+1,fp_j+1) is calculated by chain rule (see appendix).

     When the insertion of the peak birth candidate is before the first peak or after the last peak, there are no real double neighbors. In this case, we take the duplicates of peak birth candidate's unique succeeding or preceding neighbor as two virtual neighbors for proposal implementation. In this move type, we realize the sequential uniform lift for fn_*+fp_* and subsequent conditional uniform lift for fn_* within fn_* + fp_*. We claim that, the peak birth proposal by “+” move type, along with the peak death proposal by the following “–” move type, constructs a symmetric transition, i.e. equally probable events (Proposition 3). So the acceptance probability in the Metropolis-Hastings algorithm within Gibbs sampler is simply

     $$\min \left\{1,\frac{f(K,{\tilde{S}}_K,{\tilde{\sigma }}_K,\tilde{f}{n}_K,\tilde{f}{p}_K[\text{after}\hspace{0.17em}\text{peak}\hspace{0.17em}\text{death}\hspace{0.17em}\text{proposal}]|Y)}{f(K,{\tilde{S}}_K,{\tilde{\sigma }}_K,\tilde{f}{n}_K,\tilde{f}{p}_K[\text{after}\hspace{0.17em}\text{peak}\hspace{0.17em}\text{death}\hspace{0.17em}\text{proposal}]|Y))}\right\}.$$

   • 3)

     For the “–” move type, the parameter proposal process is: If K_old = K_mim we go to 1) since the lower threshold is reached; if K_old > K_mim we randomly sample one from current K_old peaks with equal probability 1/K_old, say index ^*, to delete. Then we simply abandon the associated $s_{*},{\sigma }_{*}^2,f{n}_{*}$ and $f{p}_{*}$ for likelihood reconstruction. The acceptance probability in the Metropolis-Hastings algorithm within Gibbs sampler is simply

     $$\min \left\{1,\frac{f(K,{\tilde{S}}_K,{\tilde{\sigma }}_K,\tilde{f}{n}_K,\tilde{f}{p}_K[\text{after}\hspace{0.17em}\text{peak}\hspace{0.17em}\text{death}\hspace{0.17em}\text{proposal}]|Y)}{f(K,{\tilde{S}}_K,{\tilde{\sigma }}_K,\tilde{f}{n}_K,\tilde{f}{p}_K[\text{after}\hspace{0.17em}\text{peak}\hspace{0.17em}\text{death}\hspace{0.17em}\text{proposal}]|Y))}\right\}.$$

     The “+/–” move type is demonstrated in Figure 2 for only peak location model.

     Figure 2

     Symmetric Transition for Peak Birth/death Proposal. (Each rectangle represents a uniform peak birth proposal domain, each internal vertical boundary represents a possible peak death proposal location.)

   • 4)

     For the “S” move type, the parameter proposal process is: We randomly sample one of current peaks, say ^*, with equal probability 1/K_old, and take the two neighboring peak intervals which peak ^* separates as a fused interval. A peak location is uniformly randomly drawn within this fused interval for a new candidate peak to replace s_*. The other parameters associated with this peak location mutation is kept unchanged, or they could be changed as a set within the uniformity framework. This is also a symmetric transition (Proposition 2). So the acceptance probability in the Metropolis-Hastings algorithm within Gibbs sampler is simply

     $$\min \left\{1,\frac{f(K,{\tilde{S}}_K,{\tilde{\sigma }}_K,\tilde{f}{n}_K,\tilde{f}{p}_K[\text{after}\hspace{0.17em}\text{peak}\hspace{0.17em}\text{mutation}\hspace{0.17em}\text{proposal}]|Y)}{f(K,{\tilde{S}}_K,{\tilde{\sigma }}_K,\tilde{f}{n}_K,\tilde{f}{p}_K[\text{after}\hspace{0.17em}\text{peak}\hspace{0.17em}\text{mutation}\hspace{0.17em}\text{proposal}]|Y))}\right\}.$$

   • 5)

     For the “H” move type, the posterior distribution of each ${\sigma }_k^2$ is Inverse-Gamma (ν′, η′), where $v\text{'}=v+({\sum }_{yi{j}_2\in \hspace{0.17em}peak\hspace{0.17em}K}1)/2,\text{η'}\hspace{0.17em}\text{=}\hspace{0.17em}\text{1(1/η+}{\sum }_{yij\hspace{0.17em}\in \hspace{0.17em}peak\hspace{0.17em}k}{(y_{ij}-s_k)}^2/2\text{)}$ . The joint posterior distribution of each (fn_k,fn_k, 1 – fn_k – fp_k) is Dirichlet $({\alpha }_1+n_{fn,k},{\alpha }_2+n_{fp,k},{\alpha }_3+n_{\overline{fnp{,}_k}})$ . We sample the whole set of these parameters once.

   Remark: This algorithm is a random scan (by move type probabilities) version of Gibbs sampler introduced by Gelfand and Smith (1990) with generalized parameter components:

   1. 1)

      peak (described by location, peak sample location variance, false negative and positive rates) birth or death changes the number of parameter sets by move type “+” or “–” in the aforementioned parameter sampling proposal;

   2. 2)

      peak location mutation does not change the number of parameter sets by move type “S” in the aforementioned parameter sampling proposal;

   3. 3)

      peak sample location variance or false negative/positive rate sampling does not change the number of parameter sets by move type “H” in the aforementioned parameter sampling proposal.

   The associated propositions and the detailed description of RGPMCMC are given in the appendix. The following proposition justifies the correct convergence.

   Proposition: The induced Markov chain is irreducible, aperiodic, and ergodic.

   Interpretation: In view of independent uniform proposals for non-peak parameter set and diverse move types, given any arbitrarily small neighborhood of a current state there is a positive probability that the chain lies in that neighborhood after one sampling iteration, thus the aperiodicity is verified; the irreducibility is established since the chain can move from any state to any other state one step at a time.

   3. Simulation Study

   In this section, we study the RGPMCMC performance under diverse circumstances and make comparison with other approaches in the literature.

   3.1. Prior Sensitivity Analysis for Rgpmcmc

   We consider different priors and study the discrepancy between the specified true peak locations and the estimated peak locations. Compaq Fortran 90 is our development package and we use IMSL Fortran Numerical Library to generate random numbers. 200 spectra are simulated for each of the six simulations, the set-ups and priors are listed in Table 1, where four true peaks are considered, the peak location vectors specify the true peak locations, the σ vectors specify the peak sample location variations at true peaks, the fn and fp vectors jointly specify the probabilities for us to simulate no peak sample (with probability fn), multiple peak samples (with probability fp), or single peak sample (with probability 1 – fn – fp) for each true peak from individual spectrum. The Inverse-Gamma and Dirichlet priors are for peak sample location variances and false negative and positive rates at putative true peaks. The peak number prior is K ~ U [1, 20], the starting peak number is 11, and move type probabilities are: π(+) = 0.45, π(–) = 0.45, π(H) = 0.05, and π(S) = 0.05. With burn-in 10,000 and thinning 1,000, each collection of 1,000 posterior samples takes several minutes on a PC powered by Celeron CPU. The results are given in Figure 3. In simulation 1, the true peaks are clearly clustered and the peak estimation is good. In simulation 2 with larger peak sample location variations and false positive/negative rates, the same number of peaks are identified as the true peak number under highly informative priors. By modifying peak sample location variation priors, we observe that, under large true variations, the peak estimation under less informative variation priors is worse than that under more informative variation priors. More uncertainties are introduced in simulation 3, and fewer peaks are identified than the true peak number. By modifying peak sample location variation priors, we observe that, when the variation priors are consistent with the true variation in terms of mean value, the peak number estimation seems to be better than inconsistent variation priors. Simulation 4 shows that, the informative peak sample location variation priors may lead to good estimation even when the peaks’ associated parameters are different, so does simulation 5, where unevenly distributed peaks are considered. In simulation 6, even when the parameter priors become much less informative, the peak estimation performs well, since the true peaks are sharply surrounded by corresponding peak samples. These observations show that, informative a priori knowledge is desirable for reliable estimation.

   Table 1

   Simulation Configurations (IG: Inverse-Gamma prior for σ ² , D: Dirichlet prior for (fn, fp,1 – fn – fp)).

   

   Peak location (s1, s2, s3, s4) = (1/5, 2/5, 3/5, 4/5)
   (1)	σ	0.05	0.05	0.05	0.05	(2)	σ	0.10	0.10	0.10	0.10
   	fn	0.10	0.10	0.10	0.10		fn	0.20	0.20	0.20	0.20
   	fp	0.10	0.10	0.10	0.10		fp	0.20	0.20	0.20	0.20
   		IG(5, 100), D(1, 1, 8)						IG(26, 4), D(40, 40, 120)
   (3)	σ	0.20	0.20	0.20	0.20	(4)	σ	0.02	0.08	0.08	0.02
   	fn	0.30	0.30	0.30	0.30		fn	0.10	0.20	0.30	0.40
   	fp	0.40	0.40	0.40	0.40		fp	0.10	0.20	0.30	0.40
   		IG(6, 5), D(5, 5, 5)							IG(20, 6), D(5, 5, 10)
   Peak location (s1, s2, s3, s4) = (1/16, 3/16, 7/16, 15/16)
   (5)	σ	0.02	0.08	0.08	0.02	(6)	σ	0.02	0.02	0.02	0.02
   	fn	0.10	0.20	0.30	0.40		fn	0.10	0.10	0.10	0.10
   	fp	0.10	0.20	0.30	0.40		fp	0.10	0.10	0.10	0.10
   		IG(20, 6), D(5, 5, 10)							IG(3, 2), D(2, 2, 2)

   
   Figure 3

   Simulations and Estimations (The simulation directly produces the peak samples for 200 patients which are plotted in each odd row of panels without the need for peak sample detection by data preprocessing, the y-axis is patient index and the x-axis is log(m/z). The alignments are compared by gridded walls given in each even row of panels. From top to bottom: true peaks, aligned peaks by RGPMCMC, aligned peaks by RJMCMC, aligned peaks by scale-space approach, aligned peaks by super-set approach and aligned peaks by PAM clustering algorithm.)

   3.2. Comparison with Some Non-Mcmc approaches

   We first make comparisons with other non-MCMC algorithms represented by the recently developed scale-space approach (Yu et al. 2006a), super-set approach (Yu et al. 2006b), and Partitioning around Medoids (PAM) approach (Kaufman and Rousseeuw, 1990). For the six simulations in Section 3.1, the combined peak sample locations from 200 spectra along with the optimal cluster number minimizing the “median split silhouette” (Pollard and van der Laan, 2002) are taken as PAM inputs, S-plus function pam offers the medoid (cluster center) locations. These results are also included in Figure 3. Overall, RGPMCMC performs better than non-MCMC methods; the scale-space result is little better than the super-set result in terms of robustness; the PAM algorithm is not specifically designed for peak alignment, so it may perform poorly under certain circumstances, say simulations 2 and 3, where RGPMCMC can not recover all true peaks either.

   3.3. Comparison with Reversible Jump Markov Chain Monte Carlo

   We apply the same simulated data in Section 3.1 and make comparison with reversible jump Markov chain Monte Carlo (RJMCMC) algorithm by Green (1995). RGPMCMC and RJMCMC differ in the method for proposing move type “+” and/or “–” (peak birth and/or death), where the former conditions on the active Markov chain by making use of equally probable peak birth and death proposals, while the latter makes use of additional variables to construct a one-to-one matching for dimension changing (details are given in the appendix). The same prior specifications for RGPMCMC in Section 3.1 are applied to RJMCMC. We consider two starting peak numbers, 11 and 1, both equally partitions the m/z range.

   Initial Peak Number = 11