Systematic Comparison of False-Discovery-Rate-Controlling

Article pubs.acs.org/jpr

Systematic Comparison of False-Discovery-Rate-Controlling Strategies for Proteogenomic Search Using Spike-in Experiments Honglan Li,†,∥ Jonghun Park,‡,∥ Hyunwoo Kim,§ Kyu-Baek Hwang,*,† and Eunok Paek*,‡ †

School of Computer Science and Engineering, Soongsil University, Seoul 06978, Republic of Korea Department of Computer Science, Hanyang University, Seoul 04763, Republic of Korea § Scientific Data Research Center, Korea Institute of Science and Technology Information, Daejeon 34141, Republic of Korea ‡

S Supporting Information *

ABSTRACT: Proteogenomic searches are useful for novel peptide identification from tandem mass spectra. Usually, separate and multistage approaches are adopted to accurately control the false discovery rate (FDR) for proteogenomic search. Their performance on novel peptide identification has not been thoroughly evaluated, however, mainly due to the difficulty in confirming the existence of identified novel peptides. We simulated a proteogenomic search using a controlled, spike-in proteomic data set. After confirming that the results of the simulated proteogenomic search were similar to those of a real proteogenomic search using a human cell line data set, we evaluated the performance of six FDR control methodsglobal, separate, and multistage FDR estimation, respectively, coupled to a target-decoy search and a mixture model-based methodon novel peptide identification. The multistage approach showed the highest accuracy for FDR estimation. However, global and separate FDR estimation with the mixture model-based method showed higher sensitivities than others at the same true FDR. Furthermore, the mixture modelbased method performed equally well when applied without or with a reduced set of decoy sequences. Considering different prior probabilities for novel and known protein identification, we recommend using mixture model-based methods with separate FDR estimation for sensitive and reliable identification of novel peptides from proteogenomic searches. KEYWORDS: proteogenomic search, novel peptide identification, spike-in data, simulation, false discovery rate control

■

INTRODUCTION To identify novel proteins, tandem mass (MS/MS) spectra are usually searched against proteogenomic databases, containing both reference (or known) and novel protein sequences derived from transcriptomic or genomic evidence.1 The novel protein sequences are constructed by diverse methods such as six-frame translation of whole-genome sequences,2 extracting splicing information from RNA-seq data,3 or expanding reference protein sequence databases using common or personal variant information.4 Identified novel proteins can be used to confirm or refine gene models,5 suggest novel protein-coding regions, and thus propose novel disease mechanisms 6 and open the avenue for precision and personalized medicine. Although proteogenomic searches are useful for novel protein (or peptide) discovery, the inflated nature of proteogenomic databases poses challenges in sensitive and reliable peptide identification.7 First, the inflated decoy database, accompanying the inflated target database for proteogenomic search, increases the number of high-scoring decoy hits, resulting in low sensitivity for known peptide identification.8 Second, peptides with low probabilities of © 2017 American Chemical Society

existence included in the target proteogenomic database could cause underestimation of false discovery rate (FDR) for novel peptide identification.7,9,10 To overcome these difficulties, separate and multistage FDR estimation methods were proposed for proteogenomic search.1,7,11−13 However, the accuracy of the proposed FDR estimation methods on novel peptide identification has not been assessed because it was hard to verify the existence of identified novel peptides. In other words, the sensitivity of proteogenomic search on novel peptide identification was not thoroughly evaluated. Considering that one of the main objectives of proteogenomic search is to identify novel peptides or proteins, it is of significant importance to estimate the sensitivity of novel peptide identifications. To estimate the effect of FDR controlling methods on the sensitivity and reliability of proteogenomic search for novel peptide identification, we used a public spike-in proteomic data set, obtained by spiking different amounts of 48 human proteins into a yeast cell lysate.14 We first showed that Received: January 17, 2017 Published: April 28, 2017 2231

DOI: 10.1021/acs.jproteome.7b00033 J. Proteome Res. 2017, 16, 2231−2239

Article

Journal of Proteome Research Simulated Proteogenomic Search

identifying peptides of low abundance human proteins in the spike-in experiment was very similar to identifying novel peptides from a real proteogenomic search and thus can be considered as a simulated proteogenomic search for novel peptide discovery. In this “simulated” proteogenomic search, yeast proteins were regarded as “known” and human proteins as “novel”. Among the simulated novel (i.e., human) proteins, only the 48 human proteins, which have been spiked into the yeast background, were considered true positive. In the simulated search, we evaluated three FDR estimation methods: global, separate, and multistage estimation. Each of the three methods was coupled to the target decoy search strategy15 and a mixture model-based method,16 which are two popular FDR estimators. In total, we examined six FDR control methods for proteogenomic search. Three popular database search toolsComet,17 X!Tandem,18 and MS-GF+19were tested with the six FDR control methods. First, we evaluated the degree of FDR underestimation of the six methods. Then, their sensitivities for novel peptide identification were compared. Finally, we examined the effect of decoy database size on the accuracy of the mixture model-based method. Unlike the target decoy search strategy, the mixture modelbased method can be applied without a decoy database, whose size must be the same as the target database in the target decoy method, and thus it can reduce the computational requirement by using a small-sized decoy database during search. On the basis of the experimental results, we suggest a guideline for choosing the optimal FDR control method for sensitive and reliable identification of novel peptides from a proteogenomic search.

■

We simulated a proteogenomic search by searching the MS/MS scans from the spike-in data set against a simulated proteogenomic database, containing 6752 UniProt yeast proteins (v201601), 171 121 UniProt human proteins (v201601), the 48 UPS1 human proteins, and 179 common contaminants. The sequences of the 48 UPS1 human proteins were downloaded from http://www.sigmaaldrich.com/lifescience/proteomics/mass-spectrometry/ups1-and-ups2proteomic.html. The total lengths of the protein sequences in the UniProt yeast and human (including the 48 UPS1 proteins) protein databases were 3 037 264 and 60 155 758 amino acids (AA), respectively. The yeast proteins were regarded as reference (or known) proteins. All of the human proteins were regarded as novel. Only the peptides from the 48 UPS1 human proteins were considered true positive. We converted the downloaded MS/MS files into the Mascot generic format (MGF) and the mzXML formats using the Trans-Proteomic Pipeline (TPP)21 (v.4.8.0). For simulated proteogenomic searches, three popular search toolsComet,17 X!Tandem,18 and MS-GF+19were used with the following options: 20 ppm parent mass tolerance, fully tryptic digestion, isotope error option on, and two missed cleavage sites at maximum. Carbamidomethylation at Cys was set as a fixed modification, and oxidation at Met was set as a variable modification. Additionally, 0.5 Da fragment mass tolerance was used. Real Proteogenomic Search

The MS/MS scans from the human cell line were searched against a proteogenomic database, containing 42 164 SwissProt human proteins (v201607), 528 723 three-frame translation of the human cDNA sequences from Ensembl release 84, and 179 common contaminants. The total lengths of the protein sequences in the SwissProt human protein database and the three-frame translated sequence database were 24 281 738 and 284 525 666 AA, respectively. We converted the downloaded MS/MS files to the mzXML format using TPP (v4.8.0). The converted MS/MS files were searched by Comet using the following parameters: 20 ppm parent mass tolerance, fully tryptic digestion, isotope error option on, and two missed cleavage sites at maximum. We set carbamidomethyl at Cys as a fixed modification, and methionine oxidation and N-terminal acetylation as variable modifications.

MATERIALS AND METHODS

Spike-in Proteomic Standard Data

To simulate a proteogenomic search, we used a public spike-in proteomic standard data set (PRIDE ID: PXD001819)14 downloaded from PRIDE Archive (https://www.ebi.ac.uk/ pride/archive/). The data set was obtained by spiking nine different levels of the Universal Proteomics Standard (UPS1) (Sigma-Aldrich), containing 48 human proteins, into a yeast cell lysate. We used six of the nine samples, in which 0.5, 2.5, 5.0, 12.5, 25.0, and 50.0 fmol of the 48 human proteins were spiked per microgram of the yeast background, respectively. Each sample was digested by trypsin and analyzed in triplicate by nano liquid chromatography (LC)−MS/MS using a nanoRS UHPLC system (Dionex, Amsterdam, The Netherlands) coupled to an LTQ-Orbitrap Velos mass spectrometer (Thermo Fisher Scientific, Bremen, Germany). The six MS/ MS data sets, respectively, contained 111 843, 112 777, 114 016, 119 182, 121 807, and 125 151 scans.

FDR Control

We used the target decoy search strategy (TD)15 and a mixture model-based method (MB)16 to validate database search results. In TD, an FDR estimate was calculated by dividing the number of decoy hits by the number of target hits above a score threshold. E values from the search results were used as score. To build a decoy database, reversed target-protein sequences were used. For MB, we applied PeptideProphet16 (TPP v4.8.0) for Comet and X!Tandem with the following options: semisupervised, mass accurate, minimum peptide length of 8 AA, and minimum peptide probability of 0. As for MS-GF+, the same version of PeptideProphet was not applicable. After communicating with the developer of PeptideProphet, we obtained an unpublished version of PeptideProphet. We ran the software with the parameters suggested by the developer: NONPARAM DECOY = XXX DECOYPROBS MINPROB = 0 CLEVEL = 2 ACCMASS PPM. For Comet and X!Tandem, Gaussian and Gumbel distributions were, respectively, used for modeling the

Human Cell Line Data

An MS/MS data set from a human cell line (HEK293) was searched against real proteogenomic database. The data set (PRIDE ID: PXD002395)20 was downloaded from PRIDE Archive. The MS/MS data set was obtained using an LTQOrbitrap Velos mass spectrometer (Thermo Fisher Scientific, Bremen, Germany) coupled to a high-performance liquid chromatography. The high-energy collisional dissociation method was used for peptide fragmentation, and the MS/MS spectra were obtained from the Orbitrap analyzer. The MS/MS data set contained 624 108 scans. 2232

DOI: 10.1021/acs.jproteome.7b00033 J. Proteome Res. 2017, 16, 2231−2239

Article

Journal of Proteome Research

Figure 1. Proportion of the number of novel PSMs to the number of total MS/MS spectra from the proteogenomic search results for the six spike-in standard protein mixture data sets and the human cell line (HEK293) data set. The database search was done by Comet. The PSMs of charge 2+ at 1% FDR are shown. TD: the target-decoy search strategy. MB: the mixture model-based method. GBL: global FDR estimation. SEP: separate FDR estimation. MLT: multistage FDR estimation.

a sample containing both pathogen and host proteins,22 where it is not clear to prefer one database to the other. In these cases, MLT would be not applicable.

distribution of discriminant function values for correct and incorrect peptide spectrum matches (PSMs). For MS-GF+, a kernel estimation methodthe semiparametric modewas adopted. After learning the mixture distribution, the FDR was estimated as the proportion of the expected number of incorrect matches to the total number of PSMs above a score threshold. We tested TD and MB with three FDR estimation approaches for proteogenomic search: global (GBL), separate (SEP), and multistage methods (MLT). In GBL, known and novel PSMs were used together for FDR estimation. In contrast, the FDR estimates for known and novel PSMs were separately calculated in SEP. To apply SEP with MB, a single mixture model was learned from all of the PSMs from both known and novel protein databases because it was not possible to learn a mixture model from small numbers of novel PSMs. Using the learned mixture model, the FDR estimates for known and novel PSMs were separately calculated by MB. For MLT, we performed a two-stage FDR estimation: The first stage consists of searching against the reference protein database and estimating FDR; the second stage consists of searching against the novel protein database only using scans not identified as known peptides during the first stage then applying FDR estimation. To use MLT with MB, we applied MB to the firststage search results. After the first stage, we used TD for the second-stage FDR control because we failed to learn an independent mixture model from small numbers of novel PSMs. When using MLT, the order of database search should be carefully determined. For novel peptide discovery, which is often the goal of proteogenomic search, it is natural to search against a novel sequence database after searching against a reference database first. However, there are tricky cases, such as

FDR Estimation with Varying Size Decoy Databases

Unlike TD, MB does not require that the size of decoy protein databases be the same as the size of target protein databases. We analyzed the performance of MB, with the varying size of decoy databases. Two FDR estimation methodsGBL and SEPwere used with MB for the analysis. MLT was not applicable because TD was employed as its second-stage FDR estimation (see FDR Control). Decoy databases of three different sizes were tested with MB: the original proteogenomic decoy database, the decoy database for reference proteins, and no decoy database. PeptideProphet in its unsupervised mode was employed for MB with no decoy database. Among the three database search tools, only Comet and X!Tandem were tested with varying-size decoy databases because PeptideProphet was not applicable to the search results of MS-GF+ without decoy databases.

■

RESULTS AND DISCUSSION

Comparison of Simulated and Real Proteogenomic Searches

The results of the simulated and real proteogenomic searches by Comet are summarized in Supplementary Tables 1 and 2. It is notable that the number of novel PSMs from the simulated proteogenomic searches was influenced by the level of UPS1 concentration in the spike-in data sets. The number of novel PSMs from the spike-in data set with the highest UPS1 concentration level (50.0 fmol) was 10 to 100 times larger than that from the data set with the lowest concentration (0.5 fmol). 2233

DOI: 10.1021/acs.jproteome.7b00033 J. Proteome Res. 2017, 16, 2231−2239

Article

Journal of Proteome Research

Figure 2. Relationship between the estimated and the true FDRs in (a) PSM and (b) peptide levels. Results for novel PSMs of charge 2+ from the 5.0 fmol spike-in data set are shown. The database search was done by Comet. TD: the target-decoy search strategy. MB: the mixture model-based method. GBL: global FDR estimation. SEP: separate FDR estimation. MLT: multistage FDR estimation. The black straight line is the identity line (y = x). The estimated FDRs are shown only in PSM-level because the peptide-level estimation was not applicable with MLT. For MLT, it is unclear how to filter MS/MS scans at its first stage based on peptide-level FDR filtering.

Nups > x

This trend was observed regardless of the FDR control methods used. As in the previous studies,7,12 SEP and MLT were more conservative in novel PSM identification than GBL. Interestingly, the difference in the numbers of novel PSMs between GBL and the other two methods increased as the UPS1 concentration level decreased. It suggests that the impact of FDR estimation methods on novel PSM identification is even greater when the relative quantity of novel proteins gets small. After confirming the relationship between the relative quantity of novel proteins and the number of novel PSMs from the simulated proteogenomic search, we compared the results of the simulated and the real proteogenomic searches. Because the number of PSMs in the spike-in data sets and the human cell line (HEK293) data set was not directly comparable, the proportion of novel PSMs to the number of total MS/MS spectra was compared. Figure 1 shows the comparison results for PSMs of charge 2+ at 1% FDR. Among the six spike-in data sets, the one with 2.5 fmol concentration was the most similar to the HEK293 data set in this regard. The number of novel PSMs from the HEK293 data set was affected by the FDR estimation methods similarly to the cases of spikein data sets. The number of novel PSMs identified by GBL was larger than that by SEP and MLT (see Supplementary Table 1). The results for PSMs with charge ≥3+ are summarized in Supplementary Table 2, and similar trends were observed. For the PSMs with charge ≥3+, the results from the HEK293 data set were most similar to the results from the spike-in data set with 5.0 fmol concentration, when compared by the proportion of novel PSMs to the number of total MS/MS spectra. To summarize, we observed that the spike-in data sets with UPS1 concentration levels of 2.5 and 5.0 fmol were comparable to the HEK293 data set. We examined the impact of the relative quantity of novel proteins in a sample on the sensitivity and the true FDR for novel PSM and peptide identification in the simulated proteogenomic search. The sensitivity was defined as follows

Nups where Nups is the number of PSMs (or identified peptides) from the UPS1 human proteins before FDR filtering and Nups>x is the number of PSMs (or identified peptides) from the UPS1 human proteins after FDR filtering at a score cutoff value x. The results at 1% FDR cutoff are shown in Supplementary Figures 1−4. The sensitivity of simulated proteogenomic search was largely proportional to the level of UPS1 concentration (Supplementary Figures 1 and 2). Furthermore, the sensitivities for PSMs and peptides with charge ≥3+ were zero regardless of the FDR estimation methods when the UPS1 concentration values were