MetaMSD: meta analysis for mass spectrometry data

Meta-analysis methods

MetaMSD combined differential protein results of multiple datasets (e.g., datasets from pilot and main studies) using Pearson’s test or Stouffer’s test (Fig. 1). Both tests were well-known meta-analysis tests that considered the directionality of hypotheses. Thus, it objectively resolved contradicting results between datasets. For example, a protein abundance of interest may be significantly larger in Group 1 than 2 for one dataset, but significantly smaller in Group 1 than 2 for another dataset. Based on one-sided p-value information, Pearson’s and Stouffer’s tests determined whether a protein of interest was more abundant in Group 1 or 2. Moreover, the resulting p-values/q-values reflected uncertainty of this decision. It was previously known that Pearson’s method was more sensitive to large p-values than Stouffer’s method, while Stouffer’s method was more sensitive to small p-values than Pearson’s method (Heard & Rubin-Delanchy, 2018). Researchers may choose one of these two meta-analyses approaches for their experiments depending on their study goals. Here, we elaborated both methods in terms of proteins and datasets. We let $p_{i,j}^L$ be a p-value of the one-sided (left) hypothesis test for jth protein using the ith dataset (study) where i = 1, …, K and j = 1, …, N. The number of hypothesis tests (or proteins) was denoted as N and the number of datasets (or studies) was denoted as $K\cdot p_{i,j}^L$ was a probability of falsely concluding that the jth protein in the ith dataset was less abundant in Group 1 than Group 2. Thus, as $p_{i,j}^L$ became smaller, we were more confident that the jth protein in the ith dataset was less abundant in Group 1 than Group 2. Similarly, if a right-sided p-value, $p_{i,j}^R$, was less than 0.05, then we concluded that the mean abundance of the jth protein in the ith study was significantly higher in Group 1 compared to Group 2 at the 95% confidence level. In mathematical notation, $p_{i,j}^R=1-p_{i,j}^L$. Thus, as $p_{i,j}^L$ became larger, $p_{i,j}^R$ became smaller.

Pearson’s test was a one-sided Fisher’s test. It took a minimum p-value of left- and right-sided Fisher’s tests with a Bonferroni correction. Different from Fisher’s method, it considered the directionality of the hypotheses. Specifically, Pearson’s test was based on the following procedures: (1)${\displaystyle Q_j^T=max\left(Q_j^L,Q_j^R\right),}$ where $Q_j^L=-2log\left({\prod }_{i=1}^K{p}_{i,j}^L\right)$ and $Q_j^R=-2log\left({\prod }_{i=1}^K\left(1-p_{i,j}^L\right)\right)$. And a p-value based for the jth protein was $p_j=min\left(1,2Pr\left({\chi }_{2K}^2\ge Q^T\right)\right)$. Stouffer’s test used the sum of inverse normal transformation of p-values. Stouffer’s statistics was: (2)${\displaystyle Z_j=\frac{\sum _{i=1}^K{Z}_{i,j}}{\sqrt{K}}}$ where Z_i,j was converted from a p-value from the ith dataset (or study) and the jth protein. Then, its left-sided p-value for the jth protein, $p_j^L$, was calculated based on Z_j which had a standard normal distribution under the null hypothesis. Its right-sided p-value was $p_j^R=1-p_j^L$. The final p-value was calculated by taking a side with a smaller one-sided p-value and applying a Bonferroni correction. For both Pearson’s and Stouffer’s tests, we considered all the proteins that were quantified in at least one of the studies, thus K may vary from one protein to another. If a protein was quantified in more than one study, then the meta-analysis was employed. If a protein was quantified in only one study, then the p-value from that one study was kept without employing the meta-analysis. Since we usually performed many hypothesis tests (e.g., 6,000 to 10,000 tests), we employed a false discovery procedure to correct multiple testing errors. MetaMSD reported q-values, which were the minimum false discovery rates that could be attained when calling that the abundances of proteins were significantly different between groups (Storey, 2002).

MetaMSD also measured a meta-analysis performance using the Integration-driven Discovery Rate (IDR) and the Integration-driven Revision Rate (IRR). Average IDR and IRR were reported for the simulation study. IDR represented a proportion of proteins detected by the meta-analysis that were not discovered in any of the individual studies. IRR represented a proportion of proteins detected in at least one individual study, but not by the meta-analysis. Thus, we wanted our meta-analysis method to give a higher IDR without unnecessarily increasing IRR at a given false discovery rate threshold.

MetaMSD software

MetaMSD took signs of test statistics and the corresponding p-values as input and generated meta-analysis p-values, the directionality of hypothesis tests, and the corresponding q-values as outputs (Fig. 1). Noting that major proteomic quantification software can generate test statistics and p-values, MetaMSD can be used to integrate protein quantification results generated from various proteomic software.

MetaMSD generated a tab-delimited file that contained both single dataset analyses results and a meta-analysis result. The file contained protein name, the directionalities of hypothesis tests, p-values and q-values for single analyses and meta-analysis results. It also generated several graphs/tables that described various aspects of differential proteins: (1) a plot that described numbers of detected differential proteins given q-value thresholds (Fig. 2A); (2) Summary statistics about the numbers of detected proteins and a meta-analysis diagnosis (Fig. 2B); and (3) Top-N differential proteins list detected by meta-analysis (Fig. 2C). The following command line generated a result file and graphs/tables:

The software MetaMSD, its manual, and example datasets are available at https://github.com/soyoungryu/MetaMSD under the terms of the MIT license, which allows further community driven meta-analysis development.

Simulated datasets

MetaMSD was tested using ten different simulation scenarios. For each simulation scenario, we generated one thousand simulations and calculated average performance measures of meta-analysis methods. For each simulation, multiple datasets (or studies) were generated and analyzed using MetaMSD. Each dataset contained 2n samples with n samples per group. A total of 2,000 proteins were generated per study. We varied the number of samples per group, n, where n =6, 9, or 12. Mass spectrometry was often not able to quantify all proteins in complex mixtures, and a protein quantified in one dataset might not be necessarily quantified in another dataset. Thus, we let the overlap percentage in quantified proteins between datasets to be 75% (ρ = 0.75) for the most of simulations. However, different ρ values were also explored and discussed. The data were simulated using the following schemes: (3)${\displaystyle x_{ij}\sim N\left({\gamma }_j,{\lambda }_j^2\right),}$ (4)${\displaystyle log\left({\gamma }_j\right)\sim N\left(\mu ,{\sigma }^2\right),}$ (5)${\displaystyle log\left({\lambda }_j\right)\sim N\left(\alpha +\beta log\left({\gamma }_j\right),\eta \right),}$ where x_ij was an abundance for protein j in dataset i, γ_j was a mean abundance for protein j, and λ_j was its standard deviation. In Eq. (5), we specified λ_j such that there was a linear relationship between log(λ_j) and log(γ_j). In other words, more abundant proteins had higher variance. The parameters were set as the followings: μ = 2.42, σ = 0.30, α = 1.63, β =  − 0.90, and η = 0.50. These values were obtained from label-free plasma proteomic data of mice (Kuusela et al., 2017). Raw data are available via ProteomeXchange with identifier PXD005022. The last three parameters were estimated by fitting a linear regression on log-transformed means of protein abundances and the corresponding log-transformed standard deviations. The simulated data and mouse data had similar protein abundance distributions including the mean–variance relationship. In our simulated dataset, 30% of proteins were differential proteins with 1.5, 2, or 4 fold changes in protein abundances between experiments. T-tests (with unequal variances) were performed for each dataset, then meta-analysis techniques were applied to combine resulting p-values. The following was a summary of 10 simulation scenarios:

Simulation Scenario 1: n = 6, ρ = 0.75, K = 2, α = 1.63

Simulation Scenario 2: n = 9, ρ = 0.75, K = 2, α = 1.63

Simulation Scenario 3: n = 12, ρ = 0.75, K = 2, α = 1.63

Simulation Scenario 4: n = 6, ρ = 0.75, K = 3, α = 1.63

Simulation Scenario 5: n = 6, ρ = 0.75, K = 4, α = 1.63

Simulation Scenario 6: n = 6, ρ = 0.75, K = 5, α = 1.63

Simulation Scenario 7: n = 6, ρ = 0.75, K = 2, α = 1.63 × 1.5

Simulation Scenario 8: n = 6, ρ = 0.75, K = 2, α = 1.63 × 2.0

Simulation Scenario 9: n = 6, ρ = 0.50, K = 2, α = 1.63

Simulation Scenario 10: n = 6, ρ = 0.25, K = 2, α = 1.63

Simulation Scenario 1 was used as a baseline scenario. Using Simulation Scenarios 1–3, we investigated MetaMSD performance for varying sample size (n = 6, 9, or 12). Based on Simulation Scenarios 1 and 4–6, we investigated MetaMSD performance for varying number of studies. Simulation Scenarios 1 and 7–8 were used to investigate the effects of protein expression level variability and varying experiment qualities on the meta-analysis results. Using Simulation Scenarios 1, 9 and 10, we observed how different overlap percentages in quantified proteins between studies affected meta-analysis results.

Kidney transplant proteomic datasets

Two urinary proteomic datasets of renal transplant patients were obtained from ProteomeXchange (Vizcaino et al., 2014) (PXD 002761). The details about mass spectrometry experiments can be found in Sigdel et al. (2014); Sigdel et al. (2016). In brief, these datasets contained four phenotypes of kidney transplant patients. However, in this paper, the focus was on a two-group comparison by comparing acute allograft rejection (AR) and stable allograft (STA) of kidney transplant patients. The first dataset was from an iTRAQ (Hardt et al., 2005) mass spectrometry study. Each pooled sample generated from either five patients with acute rejection (AR) or five patients with stable graft function (STA). The data were analyzed by DeconMSn (Mayampurath et al., 2008) and Sequest (Eng, McCormack & Yates, 1994). The p-values were obtained by comparing six pooled AR samples and six pooled STA samples using t-tests with unequal variances. The second dataset was from a label-free mass spectrometry study. There were 40 AR samples and 40 STA samples. However, to demonstrate the usefulness of the meta-analysis approaches in combining small- to large-scale studies, we varied the sample size for AR patients and for STA patients with n_A =6, 9, 12, or 40. The label-free dataset was analyzed by MaxQuant (Cox & Mann, 2008) and MSstats (Choi et al., 2014; Clough et al., 2012). We combined the p-values generated from the iTRAQ and label-free studies using MetaMSD.

Breast cancer proteomic datasets

Two proteomic datasets of estrogen receptor-positive (ER⁺) and triple-negative breast cancer (TNBC) breast cancer tumor samples were obtained from the following publications: Cancer Genome Atlas Network (2012) and Gámez-Pozo et al. (2015). To demonstrate the usefulness of MetaMSD in integrating medium-size datasets, we did not utilize all the patient data in these large-scale studies. Furthermore, in order to match patient’s clinical characteristics between datasets, we used the proteomic data of patients who met the following criteria: HER2-negative (Human Epidermal growth factor Receptor 2 negative) breast cancer, and Lymph node-positive with lymph node status of either N1 or N2. The first dataset (denoted as Dataset A) contained iTRAQ labeling proteomic data of 21 ER+ patients and six TNBC patients. The samples were analyzed using a nanoLC system coupled to a Q Exactive MS (Thermo Scientific, Waltham, MA, USA). The second dataset (denoted as Dataset B) contained label-free proteomic data of nine ER+ patients and nine TNBC patients. These were analyzed by a LTQ-Orbitrap Velos hybrid mass spectrometer coupled to NanoLC-Ultra system. These two clinical proteomic datasets were freely available with some restrictions. For the first dataset, protein quantification information we used in this paper, raw data, and the detailed description are freely available to the public at https://cptac-data-portal.georgetown.edu/cptac/s/S015. However, users need to agree to their data agreement policy. For the second dataset, the protein quantification information we used in this paper and detailed description about data processing are available in Gámez-Pozo et al. (2015) as Supplementary Materials. Raw mass spectrometry data are also available at Cancer Research Online (http://cancerres.aacrjournals.org/) upon the corresponding author’s approval. For both datasets, the p-values were calculated using t-tests with unequal variances and combined using MetaMSD.