Analysis of significant protein abundance from multiple reaction-monitoring data

Background Discovering reliable protein biomarkers is one of the most important issues in biomedical research. The ELISA is a traditional technique for accurate quantitation of well-known proteins. Recently, the multiple reaction-monitoring (MRM) mass spectrometry has been proposed for quantifying newly discovered protein and has become a popular alternative to ELISA. For the MRM data analysis, linear mixed modeling (LMM) has been used to analyze MRM data. MSstats is one of the most widely used tools for MRM data analysis that is based on the LMMs. However, LMMs often provide various significance results, depending on model specification. Sometimes it would be difficult to specify a correct LMM method for the analysis of MRM data. Here, we propose a new logistic regression-based method for Significance Analysis of Multiple Reaction Monitoring (LR-SAM). Results Through simulation studies, we demonstrate that LMM methods may not preserve type I error, thus yielding high false- positive errors, depending on how random effects are specified. Our simulation study also shows that the LR-SAM approach performs similarly well as LMM approaches, in most cases. However, LR-SAM performs better than the LMMs, particularly when the effects sizes of peptides from the same protein are heterogeneous. Our proposed method was applied to MRM data for identification of proteins associated with clinical responses of treatment of 115 hepatocellular carcinoma (HCC) patients with the tyrosine kinase inhibitor sorafenib. Of 124 candidate proteins, LMM approaches provided 6 results varying in significance, while LR-SAM, by contrast, yielded 18 significant results that were quite reproducibly consistent. Conclusion As exemplified by an application to HCC data set, LR-SAM more effectively identified proteins associated with clinical responses of treatment than LMM did.


Background
Discovering protein disease biomarkers is an urgently pressing issue in biomedical research [1]. Historically, the enzyme-linked immunosorbent assay (ELISA) is a highly accurate protein quantitation technique [2], representing the "gold standard" for measuring levels of specific proteins [3]. However, recent discoveries of many novel proteins, having no available antibodies, now limits the use of the ELISA method [4]. Moreover, for newly discovered proteins, development of high quality ELISA assays requires considerable time and resources [5]. Together, these shortcomings have created a need for a different technique of targeted protein quantitation [6]. One such technique, quantitative mass spectrometry of proteins has advanced significantly over the last decade, and several methods have been developed for relative quantification of targeted proteins, using this technique [7].
In recent years, multiple reaction monitoring (MRM) mass spectrometry has been developed as an attractive tool for targeted proteins and now represents a promising alternative to ELISA for quantification of proteins. To that end, MRM uses sequence-specific tandem mass spectrometer fragmentations of peptides that can provide highly selective measurements for peptide-containing proteins [8]. Thus, without enrichment or fractionation approaches, MRM assays allow for the quantitation of protein ions within ranges of low (μg/mL) to high (ng/mL) levels [5,9]. Additionally, development time for MRM assays is relatively shorter and less expensive than that for ELISA, with no requirements (and thus no costs) for antibody development.
Due to these many advantages, MRM assays are being increasingly used in systems biology and clinical investigations [10][11][12][13]. There are some analysis tools to analyze MRM data. Skyline is representative tool to create and analyze MRM dataset [14]. ProteoSign gives some simple statistical analysis results to find differentially expressed proteins [15]. roHits-viz is a web-based tool for visualizing interaction of protein data [16]. However, the development of statistical methods to determine significantly abundant proteins by the MRM assay has not received enough attention, compared to improvement of the MRM assay technology itself [17]. For MRM data analysis, two-sample t-tests or paired t-tests have been applied to identify proteins that change in abundance between two groups [18][19][20]. To test for multiple groups, one-way analysis of variance (ANOVA) has been employed [21,22].
Recently, a linear mixed model (LMM) approach was proposed for MRM data analysis, and was implemented in MSstats (v3.7.1) [17], resulting in its widespread adoption [23]. LMM approaches include a single parameter for a group effect, peptide effects, and group and peptide interaction effects, along with subject and run effects. LMMs include both subject and run effects that can be designated as either random or fixed, to identify significantly evident proteins by testing the group effect. Recently, we observed that the LMM approach often provides different p-values for the group effect, from the same data, depending on which effects are treated as random or fixed. In particularly, the LMM test results vary considerably when there is an interaction effect between group and peptide.
In this report, we propose a new logistic regressionbased method for Significance Analysis of Multiple Reaction Monitoring (LR-SAM). Unlike LMMs, our LR-SAM approach uses a much smaller number of parameters. Moreover, our LR-SAM does not require inclusion of all the effects related to the run. Accordingly, our model does not need to specify run effects as random or fixed.
Since LR-SAM uses log2-transformed relative intensity values, it does not need to include run effects. We consider two separate cases when the effects of peptides in a protein are homogeneous or heterogeneous. For the significance test for proteins, we consider Wald type tests, likelihood ratio tests (LRTs), and score tests [24].
Through various simulation studies, we compared the performance of LMMs to our approach. We first observed cases in which an LMM approach did not preserve type I error, and we also examined cases in which LR-SAM methods performed better. Some tests, based on the LR-SAM method, were more powerful when the expression pattern of peptide between groups was heterogeneous.
To further establish the translational relevance of our method, we compared LMMs to our approaches in a clinical study to identify candidate serum biomarkers of sorafenib response and prognosis in patients with hepatocellular carcinoma (HCC). Sorafenib is a multikinase inhibitor mainly used for the treatment of various solid tumors including HCC [25]. However, the response rate of sorafenib is about 8% for HCC patients. Due to its high cost, it would be more economical and more beneficial to patients if Sorafenib is applied only to the patients with high chance of response. Moreover, there are no good markers to predict the patients' response to sorafenib [26]. From May 2013 to August 2014, 115 HCC patient serum samples were collected from Seoul National University Hospital as part of an ongoing HCC study. One hundred twenty-four candidate protein biomarkers, and hepatic disease-associated proteins, were chosen from 50,265 proteins, based on the LiverAtlas database [27]. Sorafenib responses were measured according to the modified response evaluation criteria in solid tumors (mRECIST) guidelines [28]. Patients with complete responses, partial responses, or stable disease were categorized as responders, while those with progressive disease were categorized non-responders. Both LMM approaches and LR-SAM methods were employed to analyze the MRM data. peptide detected. The MRM data has a hierarchical structure such that one protein contains several peptides. Since our goal is to identify proteins that are significantly different between the group, we need to combine the summarized expression of peptides from the same protein efficiently. Also, the batch effect would occur when MRM experiments were performed in different batches separately. These characteristics of MRM make it difficult to use the standard analyses such as two sample t-test.

Review of LMM approach
For a given protein, suppose it is comprised of K peptides. Let y i, j(i), k, l denote the log2(intensity) value of the j-th subject, nested in the i-th group of the k-th peptide and the l-th run. Then, the LMM used in MSstats is given as follows: where μ is the global mean; G i is the i-th group effect; S(G) j(i) representing the j-th subject effect nested in the i-th group; P k stands for the k-th peptide effect; R l stands for the l-th run effect, (G × P) ik stands for the interaction effect between the i-th group and the k-the peptide; and (P × R) kl signifies the interaction effect between the k-th peptide and the l-th run. When all the effects are treated as fixed, these parameters have the following restrictions: P 2 i¼0 G i ¼ 0, P JðiÞ jðiÞ¼1 SðGÞ jðiÞ ¼ 0, ðP Â RÞ k;l ¼ 0; P L l¼1 ðP Â RÞ k;l ¼ 0 , and ϵ i; jðiÞ;k;l $ Nð0; σ 2 ϵ Þ. Here, G 0 stands for the effect of a reference group of MRM data. When the subject and run effects are treated as random, the restrictions on S(G) j(i) , R l , and (P × R) kl are replaced by SðGÞ jðiÞ $ Nð0; σ 2 S Þ , R l $ Nð0; σ 2 R Þ and ðP Â RÞ kl $ Nð0; σ 2 PÂR Þ, respectively. For most MRM data analyses, the investigator's interest lies in determining which proteins differ in abundance between two groups. Thus, the hypothesis of interest is given below for comparing two groups: The MSstats uses the t-test for this hypothesis.

LR-SAM with fixed effect
Consider a logistic regression model using log2-transformed relative intensity value, as follows: where Z j is a group indicator of the j-th subject that is assumed to follow a Bernoulli distribution; α is an intercept; and β k is the coefficient of the k-th peptide. Note that the β k values are related to (1). Therefore, if we treat all β k 's as fixed, the hypothesis for comparing the two groups is: To test (4), we consider the likelihood ratio test (L) with K degrees of freedom, as follows: Here,α 0 is the maximum likelihood estimate (MLE) of α under the null hypothesis (see above). l 0 ðα 0 Þ is the maximum likelihood value under the null hypothesis, andα;β 1 ; …;β K are the MLEs of α, β 1 , …, β K , respectively. It is also known that L asymptotically follows a chi-square distribution, with K degrees of freedom.
A Wald type test (W) statistic for analysis, is given below: Here, W also asymptotically follows a chi-square distribution, with K degrees of freedom under the null hypothesis. If we assume that the β k values are homogeneous, a Wald test with 1 degree of freedom (W 1 ) can be considered as follows: in whichβ p is the weighted average ofβ k as where t k is a weight given as , and A is the variance ofβ p , given as Thus, W 1 asymptotically follows a chi-square distribution, with 1 degree of freedom, under the null hypothesis. Moreover, if the β k values are homogeneous, Eq. (3) (shown above) can be reduced to the following model: The hypothesis of interest from this model (8) is H 0 : β * = 0, and a Wald test (W S ) is given by: Here,β Ã is the maximum likelihood estimate of β * from the model (8), and W S asymptotically follows a chi-square distribution, with 1 degree of freedom.

LR-SAM with random effects
When effects of peptides were heterogeneous, we could detect significant heterogenous effect by assuming random effects on coefficients β k . This assumption is commonly used in meta-analysis or rare variant analysis of genetic study [24,29]. If we assume that β k follows a normal distribution, with mean 0 and variance w k τ (where w k is a known prior weight), then the logistic regression model (3) is expanded to a mixed effect model, and the hypothesis (4) is equivalent to: Since the hypothesis H 0 : τ = 0 is on the boundary of the parameter space, the variance-component score test can be considered. The score test statistic of the variance-component for (10) is: Here,μ 0 is the estimated probability under and Y k = (y 1k , …, y nk ) ′ ; Z means the group indicator vector, and W is a diagonal matrix with the k-th element as w k .
It is known that S VC follows a mixture of chi-square distributions chi-square distributions with 1 degree of freedom, and λ k is the k-th eigenvalue of P 1/2 KP 1/2 [24]. Here, P matrix with the k-th element asμ 0k ð1−μ 0k Þ. For simplicity, we assume a flat prior weight given by w k = 1 for k = 1, …, K.

Simulation design
We next performed simulation studies to investigate the power of the LMM and LR-SAM approaches for comparing two groups, and whether they preserve type I error. For analysis, there were four LMMs available, depending on how the random or fixed effects were specified: (i) LMM(FF), with fixed subject and run effects, (ii) LMM(FR), with fixed subject effect and random run effects, (iii) LMM(RF), with random subject and fixed run effects, and (iv) LMM(RR), with random subject and run effects. For each simulated dataset, the best LMM, LMM(best), was selected among four LMMs with the smallest Akaike information criterion (AIC) value [30]. Thus, there were five LR-SAM test statistics, L, W, W 1 , W S , and S VC , for comparison. Simulation data was generated from the model, using 1000 repetitions. The global mean, μ, was arbitrarily set to 15, while the reference effects, G 0 , S(G) 0, (0) , and (G × P) 0, k , for k = 1, …, K, were set to 0. The normal distribution, with mean 0 and variance 0.5, was set as the error distribution, and the number of peptides was set to 4. For the random subject effect, we generated S(G) j(i) from the identical normal distribution independently, with mean 0 and variance 0.25 for i = 1, 2. For the fixed subject effect, we set S(G) j(i) as follows: The peptide effect, P k , was considered as a fixed effect, and was set as follows: For the random run effect, we generated R l from the identical normal distribution, independently, with mean 0 and variance 0.25. For the fixed run effect, we set R l as follows: We then considered the peptide×run interaction effect as either random or fixed, followed by the type of run effect. For the random peptide by run interaction effect, we generated (P × R) k, l from the identical normal distribution independently, with mean 0 and variance 0.1. For the fixed peptide×run interaction effect, we set (P × R) k, l as follows: Provided those parameters, we further considered four group effect scenarios (GSs), as follows: GS 0 : G i ¼ 0 for all i and GS 1 : For detecting the group×peptide interaction effect, we set (G × P) i, k , as follows: Here, using σ GP values determined from the squared average of the interaction effect, we considered three interaction effect scenarios (ISs) for σ GP , as follows:

Materials
To identify candidate serum hepatocellular carcinoma (HCC) biomarkers for prognosis and response to the tyrosine kinase inhibitor sorafenib, data from 115 HCC patients who had undergone continuous administration of sorafenib for more than 6 weeks, were collected from Seoul National University Hospital, as part of an ongoing study between May 2013 and August 2014. HCC was diagnosed by histological or radiological evaluation, with reference to the American Association for the Study of Liver Diseases (AASLD) [31] or the European Association for the Study of the Liver (EASL) [32] guidelines. All procedures/analyses were approved by the Seoul National University Hospital Institutional Review Board (IRB protocol No. 0506-150-005). Sorafenib response was evaluated using the modified response evaluation criteria in solid tumors (mRECIST) [28], using independent radiologic assessments. Patients with complete response, partial response, and stable disease were categorized as responders, while those with progressive disease were categorized as non-responders. Toward these objectives, a total of 115 serum samples were randomly separated into two batches: the 1st batch consisted of 65 samples and the 2nd batch consisted of 50 samples. Of those, the total number of responders was 40, and the number of non-responders was 75. One hundred twenty-four candidate protein biomarkers, known hepatic disease-associated proteins, were chosen from 50,265 proteins, based on the LiverAtlas database [27]. One to seven peptides, comprised within each protein, were measured, and counted (Table 2).
Since there were some cases in which LMM(FF), LMM(FR), and L did not preserve type I error in our simulation studies, we applied LMM(RF), LMM(RR), W, W 1 , W S and S VC to analyze the MRM data. Quantile normalization, provided within the MSstats package [20], was employed for preprocessing. To adjust the mean differences between the two batches, batch indicators were included in models (1), (3), and (8) for MRM data analysis.

Simulation results
The type I error and the empirical power were estimated as the proportion of p-values under 0.05, out of 1000 repetitions. The type I error rates of LMMs and LR-SAM are summarized in Table 3. The true type of subject effect more strongly affected the type I error rate of LMMs, as compared to the run effect. Analogously, when the true subject effect was fixed, the four LMMs controlled the type I errors. Among the five LR-SAM tests, L could not control the type I error rate when the sample size was 20, while the other four LR-SAM tests could. When the sample sizes were 50 and 100, all LR-SAM tests controlled type I error. When the true subject effect was random, LMM(FF) and LMM(FR) could not control type I error, whereas LMM(RF) and LMM(RR) could. L did not control the type I error rate when the sample sizes were 20 and 50, while the other four LR-SAM tests did. The AIC value of LMM(FF) tended to be the smallest among the four LMMs, under any conditions. Thus, LMM(FF) was most frequently selected as the best LMM. Some LR-SAM tests showed increased power as the interaction effects became large. When there was only an interaction effect, without group effects, the power of the LMMs W 1 , and W S did not increase, as shown in Table 4. As the sample size and the interaction effect became large, W and S VC showed increased power. Moreover, S VC provided higher power than W, when the sample size was 50, while W provided higher power than S VC , when the sample size was 100.
The size of the group effect showed large effects on all LMMs and LR-SAM methods. The estimated powers of the five LMMs, and the LR-SAM methods for GS 1, are shown in Table 5. The powers of L, W, W 1 , and S VC were all affected by the size of the interaction effect,    when the group effect scenario was GS 1, while those of the five LMMs did not. As the interaction effect increased, the powers of L, W, and S VC increased, while the powers of four of the LMMs, and W S , did not, and the power of the W 1 method even decreased. When the interaction scenario was IS 1, the powers of the L and S VC methods were higher than those of the LMMs, for a sample size of 100. For the IS 2 scenario, once the sample size exceeded 50, the L and S VC methods provided higher powers than those of the LMMs. Collectively, as the size of group effects became large, all methods provided increased power. As the interaction effect became large, only the L, W, and S VC methods provided increased power.

HCC data analysis
Using the approaches described above, the maximum -log10-transformed p-values of the linear mixed models (LMMs) were higher than those of the LR-SAM methods. However, the correlation of the transformed p-values of LMM(RF) and LMM(RR) was lower than that of LR-SAM methods, as shown in Fig. 1. The correlation between LMM(RF) and LMM(RR) was 0.6558, while the minimum correlation among the LR-SAM methods was 0.784. The correlation between W S and S VC was the highest (0.9627), among the LR-SAM methods.
The lowest p-values of LMM(RF) and LMM(RR) were lower than those of the LR-SAM methods, valued at 1.87 × 10 −18 and 1.09 × 10 −30 , respectively. The lowest p-value of W, W 1 , and W S was 2 × 10 −7 , and the lowest p-value of S VC was 3.49 × 10 −14 (Table 6). Although the lowest p-values of the LR-SAM methods were less significant than those of the LMMs, and for many proteins, LR-SAM methods provided higher significance results (Fig. 1).
The LMM(RF) and LMM(RR) methods identified 11 and 37 proteins, respectively. There were six (14.29%) proteins simultaneously identified by LMM(RF) and  (Fig. 2). All LMMs and LR-SAM methods simultaneously provided significance results for the presence of the proteins GPX3, IGHG1, IGHG3 and IGJ. Among these, IGJ (immunoglobulin J chain, linker protein for immunoglobulin alpha and mu) was previously reported as having a significant difference of expression in HCC tumors [33], while GPX3 (glutathione peroxidase 3) was reported as a tumor suppressor in HCC [34].
There were two proteins that LMMs simultaneously provided a significance result, while LR-SAM methods did not. Resultant p-values for the simultaneously identified 20 proteins are shown in Table 6. Two proteins that demonstrated significance by LMM(RF) and LMM(RR) were also previously reported: CD5L (CD5 antigen-like), reported as differentially expressed in hepatitis C patients [42]: and APOA4, (apolipoprotein A-IV), reported as misexpressed in liver metabolic disorders [43].

Discussion
We examined the performance of LMMs and LR-SAM methods, through extensive simulation studies. When the true subject effect was random, LMM(FF) and LMM(FR) did not preserve type I error (Table 3). We also applied Akaike information criterion (AIC), for model selection, to check whether or not the best-performing LMM preserved type I error. However, our empirical study showed that LMM(FF) had the smallest AIC value for most simulation settings, which made it difficult to use AIC as a model selection criterion. On the other hand, LR-SAM methods, except L, well preserved type I error under any type of subject and run effects.
Since the hypothesis (2) did not consider the interaction effect, the power of the LMM approach was not affected under any size of the interaction effect. However, when the interaction effect (without the group effect) was considered, hypothesis (2) could not perform well, as we observed through our simulation studies. On the other hand, testing hypotheses (8) with W and (10) with S VC , changes under any size group effect and interaction effect were detectable. Additionally, LR-SAM methods did not provide higher power than LMMs, when there were no or only weak interaction effects. Although LMM provided more protein significance results than LR-SAM methods, the proteins identified by the LMM(RF) and LMM(RR) approaches were very different from each other. On the other hand, LR-SAM methods provided more consistent significance results than those of the LMMs. Moreover, LR-SAM methods provided greater significance results than those identified by LMM(RF), for most proteins, and several proteins were identified only by LR-SAM methods. However, there was a still performance difference between the choices of fixed and random effect models in LR-SAM methods, as was in the LMM. This difference was caused by different testing hypotheses between W and S VC in LR-SAM. In addition, note that S VC allows the heterogenous effects of peptides, while W does not.

Conclusion
LMMs have been widely used to identify proteins with significantly altered abundance, in distinct disease states, based on MRM assays [39]. However, we found that LMM approaches provide inconsistent significance results for the same MRM data, depending on which effects are treated as random or fixed by simulation results. It is a well-known property of LMMs that the variance of model parameters are underestimated, when the fixed effect model is fitted but the true effects are random, and vice versa [44]. As a result, the protein significance results of LMMs may vary, depending on whether the true effect is random or fixed, and we also observed this phenomenon, as shown in Fig. 1 and Table  6. Thus, it is highly important to correctly specify the effect as random or fixed.
In this paper, we propose a new logistic regression-based method for Significance Analysis of Multiple Reaction Monitoring (LR-SAM). Unlike LMMs, our LR-SAM approach uses a much smaller number of parameters. Moreover, our LR-SAM does not require inclusion of all the effects related to the run. Accordingly, our model does not need to specify run effects as random or fixed.
In simulation study, our LR-SAM preserved type I error when the true subject effect was random, while LMM(FF) and LMM(FR) did not. In real data analysis, although LMM provided more protein significance results than LR-SAM methods, LR-SAM methods provided more consistent significance results than those of the LMMs. Thus, our proposed method, LR-SAM could give more reliable results than previous protein studies.