Comparison of False Discovery Rate Control ... - ACS Publications

Mar 20, 2017 - approach resulted in an increased number of false identifications. ... proteomics, variant peptides, false discovery rate, cancer prote...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/jpr

Comparison of False Discovery Rate Control Strategies for Variant Peptide Identifications in Shotgun Proteogenomics Mark V. Ivanov,†,‡ Anna A. Lobas,†,‡ Dmitry S. Karpov,§,∥ Sergei A. Moshkovskii,§,⊥ and Mikhail V. Gorshkov*,†,‡ †

Institute for Energy Problems of Chemical Physics, Russian Academy of Sciences, Moscow 119334, Russia Moscow Institute of Physics and Technology (State University), Moscow Region, Dolgoprudny 141700, Russia § Institute of Biomedical Chemistry, Moscow 119121, Russia ∥ Engelhardt Institute of Molecular Biology, Russian Academy of Sciences, Moscow 119991, Russia ⊥ Pirogov Russian National Research Medical University, Moscow 117997, Russia ‡

S Supporting Information *

ABSTRACT: Proteogenomic studies aiming at identification of variant peptides using customized database searches of mass spectrometry data are facing a dilemma of selecting the most efficient database search strategy: A choice has to be made between using combined or sequential searches against reference (wild-type) and mutant protein databases or directly against the mutant database without the wild-type one. Here we called these approaches “all-together”, “one-by-one”, and “direct”, respectively. We share the results of the comparison of these search strategies obtained for large data sets of publicly available proteogenomic data. On the basis of the results of this evaluation, we found that the “alltogether” strategy provided, in general, more variant peptide identifications compared with the “one-by-one” approach, while showing similar performance for some specific cases. To validate further the results of this study, we performed a control comparison of the strategies in question using publicly available data for a mixture of the annotated human protein standard UPS1 and E. coli. For these data, both “all-together” and “one-by-one” approaches showed similar sensitivity and specificity of the searches, while the “direct” approach resulted in an increased number of false identifications. KEYWORDS: proteogenomics, shotgun proteomics, variant peptides, false discovery rate, cancer proteome



INTRODUCTION Recent advances in the technologies for deep proteome analysis using high-resolution mass spectrometry have stimulated comprehensive characterization of cancer cells at the protein level.1,2 In the related “proteogenomics” approach,3 one of the main objectives is identification of proteins with point mutations using customized databases consisting of the reference (wild-type) and mutant forms of proteins obtained from either DNA or RNA sequencing. The approach involves a computational workflow for identification of variant peptides originating from the mutant parts of the proteins. This approach becomes a significant part of cancer research focused on revealing the translation of somatic mutations into protein expression as more genomic data are becoming publicly available.4−6 One of the problems associated with the identification of variant peptides is the increasing size of the mutant protein database. This results in increasing number of false variant identifications when both wild-type and variant peptides are being filtered together to a desired level of false discovery rate (FDR). This further raises a question of correct estimation of © 2017 American Chemical Society

the FDR for variant peptides, a problem acknowledged in a number of recent publications.7−12 Indeed, while FDR for all identifications (“global” FDR) is close to the one estimated using target-decoy approach (TDA), it can be strongly biased to higher values for the selected group of variant peptides.13,14 The problems with FDR estimation and the identification efficiency for variant peptides in proteogenomic studies are closely related to the choice of the database search strategy. In general, there are three strategies used for the estimation of FDR for variant peptides. In the first one called here the “alltogether” approach and also known as “Separate-FDR”,7 MS/ MS spectra are searched against a combined database consisting of both wild-type and mutant proteins. In the second strategy called here the “one-by-one” approach, or “two-stage-FDR”,7 MS/MS spectra are first searched against the wild-type database, followed by the search of nonsignificant spectra against the variant one. Finally, in the third one, called here the “direct” approach, MS/MS spectra are searched against a Received: November 29, 2016 Published: March 20, 2017 1936

DOI: 10.1021/acs.jproteome.6b01014 J. Proteome Res. 2017, 16, 1936−1943

Article

Journal of Proteome Research database consisting only of mutant proteins.10 All strategies may provide different numbers of false identifications and may have different efficiencies. The difference in efficiencies may become crucial for the outcome of proteogenomic studies because the number of typically identified variant peptides can be very small ranging from several peptides to ∼100 peptides compared with tens of thousands of wild-type peptide identifications from the same data.7,11−13 While the “one-byone” approach may be advocated by the biologically sound assumption that mutant proteins and their wild-type counterparts should have different probabilities of being present in a sample, this approach may result in the loss of many if not all variants. They may simply be buried among the false-positive identifications in the first round of search against the wild-type database. Indeed, needless to say, the number of false-positive identifications of wild-type peptides may easily exceed the total number of identified variants by the order of magnitude for a typical deep proteogenomic analysis of a cancer cell line. On the contrary, the “direct” approach seems free from this drawback and sounds like the most efficient way to identify as many variants as possible for a given data set. However, the obvious danger of this approach is an artificial increase in the probability that a particular MS/MS spectrum originating from a wild-type peptide or its modified form will match the wrong variant sequence candidate due to the lack of alternatives. In this study we evaluated all three strategies to answer the question of proper selection of the variant search strategy for the proteogenomic analysis. We compared the above approaches using publicly available LC−MS/MS data sets from a number of recent large-scale proteogenomic studies. Note that while the advantages of the “one-by-one” approach have been discussed recently by Woo et al.,12 the reports that involve data-extensive comparison of all three approaches are still missing in the literature. Importantly, the results of this evaluation may be extended beyond the proteogenomic area and can be applicable in a number of proteomic or metaproteomic applications when a small subset of proteins distinct from the core sample proteome needs to be identified, such as the presence of viral proteins, proteins with specific post-translational modification (PTMs), and so on.



elsewhere5 and downloaded from CellMiner (http://discover. nci.nih.gov/cellminer). Finally, to validate the results of comparison of the search strategies, we extended our analysis using LC−MS/MS data obtained for the mixture of annotated universal protein standard UPS1, containing 48 recombinant human proteins (http://www.sigmaaldrich.com/life-science/proteomics/massspectrometry/ups1-and-ups2-proteomic.html), and E. coli proteins.15 In the cited work the tryptic peptides of UPS1 proteins were spiked into the E. coli tryptic digest, followed by deep proteome analysis. Two different protein databases, human and E. coli (described in more detail below), were used to evaluate the search strategies. This latter analysis presents a reconstruction of the variant peptide search in a typical proteogenomic study. Because the UPS peptides are known a priori, this experiment can be considered as a control for the results obtained here for the large-scale proteogenomic data. Search Parameters

Two search engines, X!Tandem (ver. CYCLONE 2012.10.01.1)16 and MSGF+ (v. 9079)17 were used to analyze LC−MS/MS data. A Python tool pepxmltk (available at https://pypi.python.org/pypi/pyteomics.pepxmltk) was used to convert X!Tandem output files to open pepXML format.18 For all searches, carbamidomethylation of cysteine and methionine oxidation were set as fixed and variable modifications, respectively. Enzyme specificity was set to “trypsin”. For the colon cancer data two missed cleavages were allowed, precursor mass accuracy was set at 10 ppm, and fragment ion mass measurement accuracy was 0.3 Da. For the NCI-60 data, two missed cleavages were allowed, precursor mass accuracy was set to 15 ppm, and fragment mass tolerance was 0.03 Da. For the UPS/E. coli data, one missed cleavage was allowed, precursor mass accuracy was set in the range of −10 to +15 ppm (due to the varying systematic error), and the fragment mass tolerance was 0.015 Da. The search parameters for the MSGF+ were the same as the ones used for colon cancer data, except two parameter options. MSGF+ does not have a limitation for maximal number of missed cleavages and also does not have a numerical value for fragment mass tolerance, and thus the value “low-res LCQ/LTQ” was used instead of 0.3 Da for the fragment mass tolerance.

EXPERIMENTAL SECTION

Sequence Databases

Data Sets

To explore the impact of the database construction on results of data set analysis, two typical approaches for customized mutant protein database generation have been considered in this study: (1) The customized variant database was generated for the analyzed sample and contained sample-specific mutations only and (2) the variant database was extended by inclusion of mutated proteins from the other samples. The first approach was used for the colon cancer data, for which the variant database was split into smaller ones, each one being specific for a given patient sample. The second approach was considered for the analysis of data from nine NCI-60 cell lines, in which the variant database contained mutant proteins from all 60 cell lines of the panel. Both mentioned approaches for database construction were used for the search of UPS/E. coli data. In particular, the UPS+E. coli database was representing the analogue of sample-specific variant database, and the Human+E. coli database was constructed as an analogue of extended variant database containing non-sample-specific variants.

We used publicly available MS/MS data containing over 13 million high-resolution spectra obtained using LTQ Orbitrap mass spectrometer for 95 samples from 90 patients with colorectal cancer.1 Because RNA-seq data on point mutations and, consequently, the customized protein databases are available for 86 samples, only those were further processed for the comparison of search strategies. The customized RNAseq-derived protein databases for each of the 86 samples contain from 3540 to over 10 000 protein sequences. The other publicly available data sets considered herein have been obtained in a large-scale deep proteome characterization of NCI-60 cancer cell lines.2 The data used in our evaluation include 2.5 million high-resolution MS/MS spectra of peptides from nine cell lines and were obtained using 1D SDSelectrophoresis prefractionation of the protein extracts, followed by in-gel digestion and high-resolution LC−MS/MS analysis on Thermo Orbitrap Elite mass spectrometer. The exome sequencing data for these nine cell lines were obtained 1937

DOI: 10.1021/acs.jproteome.6b01014 J. Proteome Res. 2017, 16, 1936−1943

Article

Journal of Proteome Research FDR Filtering

The customized FASTA databases generated using mutations predicted from RNA-seq1 for the colon cancer tissue samples were divided into wild-type (33 731 protein sequences) and variant databases (3540 to 10 018 sample-specific protein sequences). Decoy databases containing reversed protein sequences were used in all searches. Note that although all reversed, shuffled, or randomized approaches for decoy database generation are known to provide similar efficiencies,19 employing shuffled or random protein sequences is not applicable for the TDA in the case of mutant proteins. Indeed, the variant proteins share a major part of the sequence with respective wild-type proteins; therefore, adding a mutant protein to a database increases the set of target peptides in few instances only. This is also the case for decoy peptides if those are generated by reversing the target sequences. On the contrary, shuffling or randomizing the sequences of decoy proteins produces unique decoy peptides for all target proteins, even if they have substantially shared sequences. A simple example of the mentioned difference between reversed and shuffled mutant proteins is shown in Figure 1.

All nonvariant peptide-spectrum matches (PSMs) were excluded from the analysis before the FDR filtering. PSMs were considered variant-related if all protein candidates came from the variant database. Thus the peptides from nonmutated parts of the protein sequences having both mutant and wildtype proteins in the protein lists were excluded. Note that custom exclusion of wild-type peptides is required not only for the “all-together” approach but also for the “one-by-one” and “direct” approaches. The main reasoning here is that mutant databases usually contain the full protein sequences instead of the variant peptides, and thus the wild-type peptides from nonmutated parts of the proteins are still presented in the results, which is discussed in detail below. The FDR estimations for the variant peptide identifications were based on targetdecoy method (TDA),21 and the following equation for FDR was used19 FDR =

number of variant decoy PSMs number of variant target PSMs

(1)

The FDR level of 1% was set at the PSM level, and the identifications were postsearch validated using MP Score22 algorithm for X!Tandem and Percolator23 (v3.01) for MSGF+. In the “one-by-one” approach applied to the UPS/E. coli data, the first search was done against the E. coli database, followed by the search against the “customized” UPS database (Table S1) after removing all MS/MS spectra identified at 1% FDR in the first search. To further extend the comparison of the strategies we also used the whole human protein database in the second search. This resembles the situation typically occurring in the proteogenomic analysis when the variant database significantly exceeds the sample-specific variant sequence database generated using exome or RNA-seq data. For example, this would be the case when a search for the variant peptides in a probe of a patient is performed against the variant database generated using exomes of a larger group of patients with the same type of cancer.

Figure 1. Scheme of tryptic peptides generation for reversed and shuffled mutant proteins. While reversed decoy proteins produce the same number of tryptic peptides as the target proteins, the shuffled proteins produce a higher number of decoy tryptic peptides.



RESULTS AND DISCUSSION

Hidden Pitfalls of Proteogenomic Approaches to Variant Peptide Identification

For the NCI-60 data set the reference database UniProt Human (ver. 2013_09, contains 88 277 protein sequences, downloaded from http://www.uniprot.org/taxonomy/ complete-proteomes) and the NCI-60 variant database previously described14 were used. Human (20 193) and E. coli (4431) proteins from the SwissProt database were used for the UPS/E. coli data set. The UPS proteins (48 sequences) and reliably identified contaminants (10 sequences, mostly keratins) were assembled into a separate database. The list of the proteins in this database is shown in Table S-1. For the analysis of the wild-type database size influence, all proteins from Human, E. coli, and those sharing at least one tryptic peptide of length between 5 and 50 amino acid residues with any human protein were excluded from the full SwissProt database. The total number of proteins in the resulting database was 276 883. E. coli protein database was extended by copies of the generated database containing 10, 50, and 100% of these 276 883 sequences and used for the analysis as wild-type database. All data mining and manipulations were performed using inhouse developed Python scripts based on the Pyteomics library.20

The proteogenomic data analysis should be done differently from the one routinely employed in proteomics. Indeed, it has unique features and hidden pitfalls that should be well understood. The first potential pitfall may appear during the search against combined database (reference + mutant), followed by the filtering of the identifications from both databases together.7 Consider the results from combined database search having 100 000 PSMs at 1% FDR, 90% (90 000) coming from the first database and 10% (10 000) from the second one. At the same time, assume also that the second database is nine times larger. The estimated number of false-positives for the whole set of identifications is 1000. Now, because the false identifications are random, we should expect that among 1000 false identifications, 900 originate from the second database and the remaining 100 originate from the first one. Thus despite the global FDR of 1% the FDR for the identifications from the first and the second databases will be 0.11 and 9%, respectively. The above case is typical for the searches of variant peptides using the mutant protein databases, the size of which may be significantly larger than the size of the reference (wild-type 1938

DOI: 10.1021/acs.jproteome.6b01014 J. Proteome Res. 2017, 16, 1936−1943

Article

Journal of Proteome Research protein) database, yet the number of identified variant peptides is significantly smaller compared with the total number of identifications. As a result, the real FDR for the subset of variant peptides becomes significantly >1%. Thus the correct way to obtain the results with desired FDR ratio for the variant peptides using TDA is to filter only the variant PSMs, excluding all wild-type PSMs, as reflected in eq 1. The second problem similar to the one described above may occur when the posterior error probabilities (PEPs) are used for the FDR estimation. Typically, PEPs are calculated using the functional approximation of distributions of false and true identifications. However, this approach is not working when one considers a specific group of identifications due to the differences in distributions of wild and variant identifications, as shown in Figure S-1. The same scores may result in different PEP values depending on the group of chosen peptides. Thus we do not suggest employing the PEP values generated using both wild-type and variant identifications together. Finally, special attention should be paid to a search employing “one-by-one” approach when the databases share some peptides. If this is the case, the true identifications from the second search will most probably be unique for the second database. Indeed, in the “one-by-one” strategy, the spectra reliably matching the sequences from the first database are removed after the first search. The spectra matching false identifications will contain sequences from both databases. If the sequences from the first database are not removed before the filtering of the results of second (variant) search, then the FDR for the subset of identifications from the second search will be significantly overestimated. An example of this situation is shown in Figure S-2 for the search results obtained in this work. As shown in these Figures, the FDR threshold becomes too conservative when the shared peptides are not excluded before the filtering. This results in a decrease in the number of variant peptide identifications.

Figure 2. Workflow for the “one-by-one” search strategy employed in this work. Validated identifications from the first search against the wild-type database were excluded from the second search against the mutant protein database.

number of peptides in the search space of the “all-together” approach is eight times higher compared with the second search using the “one-by-one” approach. Also, the e values reported by X!Tandem search engine are used for postsearch validation. However, X!Tandem calculates and also reports raw search score called hyperscore. The hyperscore depends on the intensities and number of matched fragment ions and is independent of the search space. Therefore, ranking identifications using either e values or the hyperscores may produce different effects on the results of the “one-by-one” strategy. To find out which of the sources of variant peptide identification deficiency has greater impact, we processed the colon cancer data set using both e values and hyperscores. Figure 3 shows the results of comparison of the two search strategies. The ratio between the numbers of variant peptides identified by “all-together” and “one-by-one” approaches exceeded 1 for all data sets. This unambiguously demonstrates the advantage of the former approach. On average, it gave 11, 8, 5, and 4% more variant PSMs at 1% FDR for the colon cancer, NCI-60, human/E. coli, and UPS/E. coli data, respectively. Importantly, the results obtained for the UPS/E. coli data set used for two approaches for database construction have shown that the difference between “all-together” and “one-by-one” strategies is independent of the size of “variant” database. Expectedly, the search with smaller database containing only UPS proteins has resulted in the increase in the absolute number of identified PSMs from UPS proteins. The detailed results for all data sets are summarized further in Table S-2 and Table S-3. To reveal the possible origin of the lower efficiency of the “one-by-one” strategy, the additional variant PSMs identified in “all-together” search and missed in the “one-by-one” method were analyzed in more detail. The spectra of these PSMs may be identified as the wild-type peptides in the first stage of “oneby-one” search. Hence, they do not pass into the second stage of the search against the variant database. Alternatively, these spectra match the same variant peptides in the second search, but they get a different (typically lower) score due to the changes in the search space and thus cannot pass the FDR

“One-by-One” versus “All-Together” Approaches

The two approaches to variant peptide identification were evaluated using three data sets described in the Experimental Section: 86 TCGA colon cancer samples, 9 cell lines from NCI60 panel, and UPS/E. coli mix. The “all-together” approach implied a single search against a merged wild-type and customized sequence database. The “one-by-one” approach is schematically illustrated in Figure 2. As discussed above, the workflow shown in this Figure depicts a possible drawback of this approach. The MS/MS spectra giving false-positive identifications in the first search against the reference database are excluded from the second search against the customized database. However, these spectra may also correspond to variant peptides. The absolute number of spectra matched to false-positives may be significant compared with the total number of identified variant peptides, even at FDR of 1%. On the contrary, the “all-together” approach is free from this drawback, and this may result in its higher performance in identification of variant peptides. In other words, while the “alltogether” strategy implies, to some extent erroneously, equal a priori probabilities for variant and reference peptides to match their respective MS/MS spectra during the search, the “one-byone” strategy simply puts zero probability for many of the variant peptides to be presented in the sample. Typically, the search engines report statistical scores, such as e values or p values depending on the search space used for identification. For example, in the case of colon cancer data, the 1939

DOI: 10.1021/acs.jproteome.6b01014 J. Proteome Res. 2017, 16, 1936−1943

Article

Journal of Proteome Research

Figure 3. Results of applying “all-together” and “one-by-one” strategies for different data sets: (a) colon cancer, (b) NCI60, and (c) UPS/E. coli with UPS database used as a “customized” database mimicking mutated proteins and (d) UPS/E. coli with whole human database used as the “mutant” database. The efficiencies of the strategies are compared by the numbers of identified variant PSMs (or UPS PSMs for the UPS/E. coli data set).

only due to the false identifications originating from the first search in the “one-by-one” strategy. When colon cancer identifications were filtered using hyperscore instead of e value, the fraction of variant PSMs shared between the two approaches reached 100% for most samples. Because X!Tandem’s hyperscore does not depend on the search space, the variants receive the same scores for both “alltogether” and “one-by-one” strategies. At the same time, the advantage of the “all-together” approach decreased from 11 to modest 3% (see Table S-2, Table S-3, and Figure S-3).

threshold. The results of this analysis are shown in Figure 4. Indeed, we found that by median, 80, 75, 40, and 10% of variant

Recall-Precision Plot for UPS/E. coli Data Set

Using one of the UPS/E. coli experiments, we estimated sensitivity and specificity of three approaches under study and plotted the recall-precision plot. All E. coli peptides were filtered out from the results. Recall was calculated as the fraction of identified UPS peptides to all UPS peptides in the database. Precision was calculated as the fraction of identified UPS peptides to all identified human peptides, assuming all human peptides not related to UPS proteins or contaminants being false matches. The resulting recall-precision plot is shown in Figure 5a; “all-together” and “one-by-one” approaches show similar sensitivity and specificity with the area under recallprecision curve (AUC) equal to 0.902 and 0.895, respectively. However, the search strategy called here “direct” shows the worst result with AUC of 0.843. Thus the use of “direct” search strategy leads to increased number of false-positive identifications. That can be explained by the homology of human proteins and E. coli proteins excluded from the search, as discussed by Noble et al. in the already mentioned work.10 This

Figure 4. Analysis of the extra PSMs identified using the “all-together” strategy. A significant fraction (shown in the box plots) of these variant PSMs matched wild-type peptides in the first stage of the “one-by-one” search.

PSMs from the UPS/E. coli, human/E. coli, NCI-60, and colon cancer data, respectively, that were missed by the “one-by-one” search and identified using the “all-together” strategy, match the wild-type sequences during the first round of “one-by-one” search. Note also that the fractions of variant PSMs shared between “all-together” and “one-by-one” results are similar for all data sets but smaller than 100% intersection. The latter would be expected if the difference between the strategies was 1940

DOI: 10.1021/acs.jproteome.6b01014 J. Proteome Res. 2017, 16, 1936−1943

Article

Journal of Proteome Research

the variant peptides from the second search was kept constant at 1% level. As shown in Figure 6, increasing FDR threshold

Figure 5. Recall-precision (a) and q-value (b) plots obtained using “all-together”, “one-by-one”, and “direct” strategies for UPS/E. coli data set. For the recall-precision estimation the UPS peptides and other human peptides were considered as relevant and nonrelevant elements, respectively. Only human PSMs were accounted for in the qvalue estimation.

Figure 6. Number of identified variant PSMs at 1% FDR level for different search strategies. Color bars represent different approaches evaluated in this study: green, “all-together”; red, “direct”; blue, “oneby-one”. The latter approach was evaluated for a range of FDR levels set for the search against wild-type database.

homology leads to high-score matches of human non-UPS peptides with sequences close to E. coli peptides really present in the sample. Going further in discussion of the results of comparison, note that the “direct” approach was also left behind in terms of the total number of human PSMs, as shown by q-value plot in Figure 5b. The “direct” approach results in fewer identifications when FDR is set in a range of 0 to 2%. For the FDR of 2% and higher it starts getting close in efficiency to the “one-by-one” approach. This happens because the spectra from E. coli peptides are matched by false identifications in the “direct” approach, while in the “all-together” and “one-by-one” approaches these spectra are matched by true E. coli peptides. As an example, consider the decoy identification with highest score in the “direct” search, which is LAGLGGDEPGGRGALR peptide from human decoy database. However, the underlying MS/MS spectrum is matched with better score by E. coli target peptide ALGLNDELAHSSIR in the case of both “all-together” and “one-by-one” approaches. Simply speaking, a researcher performing blind data analysis does not know which variant peptides are true or false. In this case, the “direct” approach may report the same number of identifications as the other strategies at the FDR level of 2% and higher. However, because the search database did not contain the peptides really present in the sample, false variants received their chance to match more PSMs. This is a typical example of the fundamental problem with the search strategies employing only targeted databases, which do not contain all possible sequences present in the sample. Furthermore, in this work, due to the persistently lower number of identifications and the problems described above, the “direct” search strategy was not used for the proteogenomic data sets.

resulted in significant decrease in the total number of variant identifications. One could expect that 0% FDR used for the first search would provide the best results for the second search against the variant database; however, in this case one will face the problem similar to the one discussed above for “direct” searches. The wild-type target false or low-score true identifications are replaced with variant decoy peptides, which increases the threshold for variant PSMs filtering (see Table 1). Table 1. X!Tandem Identifications of Spectrum A0218_4I_R_FR14.2512.2512.2 in TCGA Data Set for Different Search Strategies

sequence hyperscore e value decoy? variant?

all-together

wild-type

one-by-one, 0% wild-type FDR

one-by-one, 1% wild-type FDR

HLPTIR 24 0.16 no no

HLPTIR 24 0.027 no no

HLTLPR 19.3 0.0096 yes yes

unmatched − − − −

While nonvariant target peptide HLPTIR matches the corresponding MS/MS spectrum better than decoy variant peptide HLTLPR, the HLPTIR identification cannot pass 0% FDR threshold in the first search against wild-type database, and the spectrum is passing to the second search, where it is matched by decoy variant. This example speaks further against the use of 0% FDR in the multistep search strategies. Effect of the Size of the Wild-Type Protein Search Space

The size of the wild-type protein database may be a factor affecting the efficiency of one or the other search strategy. To evaluate this effect, the databases with different numbers of SwissProt proteins from different organisms were added to E. coli database, as described in the Experimental Section. Here we generated four “wild-type” databases containing E. coli (4431 proteins), E. coli + 10% of SwissProt (32 119), E. coli + 50% of SwissProt (142 872), and E. coli + 100% of SwissProt (281 314), respectively. The q-value and recall-precision plots obtained for the “one-by-one” approach are very similar for

Closer Look at the “One-by-One” Strategy: Selecting the FDR Thresholds for the First Search

The idea that the “one-by-one” approach reports fewer identifications than “all-together” due to loss of variants among the false identifications in the first search against the wild-type database was probed by varying the FDR filtering level for the first search. In this evaluation we considered one of the colon cancer data sets, for which the FDR of the first search against wild-type database was varied up to 50%. The FDR for 1941

DOI: 10.1021/acs.jproteome.6b01014 J. Proteome Res. 2017, 16, 1936−1943

Article

Journal of Proteome Research

Figure 7. Total number of variant PSMs at 1% FDR for different number of PSMs used for Percolator training. Green and blue bars show results for the “all-together” approach, and the red bar shows results for the “one-by-one” approach.

all four (see Figure S-4). For the “all-together” approach the qvalue plots are also similar; however, the recall-precision plots show a small advantage of using larger size wild-type database (Figure S-5). That behavior can be explained by the filtering procedure used here, when we filter out all nonhuman identifications (analogue of the wild-type PSMs) and leave only the human ones (analogue of the variant PSMs), and because of increasing nonhuman part of database, the false human PSMs are replaced with false PSMs belonging to other organisms.

is obviously the fastest because the search database is decreased by wild-type proteins.



CONCLUSIONS Our comparison of three possible variant peptide search strategies used in proteogenomics has shown that the “direct” approach based on the search against target protein sequences only carries a danger of providing a large number of false variant identifications. Moreover, there is no feasible way to control the presence of these false identifications in the search results. At the same time, the “all-together” and “one-by-one” strategies demonstrate similar performances in terms of sensitivity and specificity with the former one providing more identifications for all data sets evaluated. However, if the workflow is based on a search engine relying on a large size search space for correct score calculation, for example, X!Tandem and its expectation value score, the “alltogether” strategy should be used. The same is true if the workflow also includes a postsearch validation based on the data training, such as Percolator. The results obtained in this study advocate for using the “alltogether” search strategy in proteogenomics. However, similar efficiency can be achieved using “one-by-one” (or “two-step”) strategy if the potential caveats of this approach such as the effect of search database size or the size of the training data set used for postsearch validation are properly taken into account.

Switching the Search Engine and the Validation Algorithm: Evaluation of the Search Strategies Using MSGF+ and Percolator

The use of search engines employing principally different scoring algorithms and postsearch validation tools may also affect the results of the variant search using different strategies. Data sets of colon cancer were additionally processed using MSGF+ search engine and Percolator. We found that the “alltogether” approach produced on the average 11% more variant PSMs at 1% variant FDR (see Table S-4). One of the possible explanations for such a gain is that Percolator uses true identifications for training underlying SVM algorithm. The “alltogether” strategy provides a larger number of identifications for training because of the presence of wild-type peptides. To support this assumption, we varied a number of PSMs used by Percolator for training by changing its input parameter called “--subset-max-train” (see Figure 7). The decrease in the number of training PSMs to 500 resulted in closer efficiencies of “all-together” and “one-by-one” approaches.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jproteome.6b01014. Supplementary Figure 1S. PEP calculation for an example of combined database search. Supplementary Figure 2S. Distributions for the X!Tandem e values for the PSM identifications obtained for TCGA-A6-380701A-22 sample from the colon cancer data set using “one-by-one” strategy. Supplementary Figure 3S. Percentage of variant identifications shared between the

Computational Efficiency

There is no significant difference in the computational efficiencies between “all-together” and “one-by-one” approaches. For example, the analysis of a single colon cancer data file using MSGF+ and Percolator took 19 and 24 min for “all-together” and “one-by-one” strategies, respectively. However, the latter approach requires more time from a researcher to manage the procedure of removing the spectra identified in the first stage of wild-type protein search. The “direct” strategy 1942

DOI: 10.1021/acs.jproteome.6b01014 J. Proteome Res. 2017, 16, 1936−1943

Article

Journal of Proteome Research



results obtained using “all-together” and “one-by-one” strategies for all data sets filtered hyperscore. Supplementary Figure 4S. Recall-precision and q-value plots calculated using “all-together” strategy for UPS1 data set with different percentages of added SwissProt proteins to the E. coli (“wild-type”) database. (PDF) Table 1S. List of UPS proteins and reliably identified contaminant proteins. (XLS) Table 2S. Number of identified variant PSMs at 1% variant FDR by X!Tandem for all used data sets. (XLS) Table 3S. Number of identified variant peptides at 1% variant FDR by X!Tandem for all used data sets. (XLS) Table 4S. Number of identified variant PSMs at 1% variant FDR by MSGF+ and Percolator for colon cancer data. (XLS)

(8) Nesvizhskii, A. I. Proteogenomics: concepts, applications and computational strategies. Nat. Methods 2014, 11 (11), 1114−1125. (9) Chernobrovkin, A. L.; Zubarev, R. A. Detection of Viral Proteins in Human Cells Lines by Xeno-Proteomics: Elimination of the Last Valid Excuse for Not Testing Every Cellular Proteome Dataset for Viral Proteins. PLoS One 2014, 9 (3), e91433. (10) Noble, W. S. Mass spectrometrists should search only for peptides they care about. Nat. Methods 2015, 12 (7), 605−608. (11) Lobas, A. A.; Karpov, D. S.; Kopylov, A. T.; Solovyeva, E. M.; Ivanov, M. V.; Ilina, I. Y.; Lazarev, V. N.; Kuznetsova, K. G.; Ilgisonis, E. V.; Zgoda, V. G.; et al. Exome-based proteogenomics of HEK-293 human cell line: Coding genomic variants identified at the level of shotgun proteome. Proteomics 2016, 16 (14), 1980−1991. (12) Woo, S.; Cha, S. W.; Bonissone, S.; Na, S.; Tabb, D. L.; Pevzner, P. A.; Bafna, V. Advanced proteogenomic analysis reveals multiple peptide mutations and complex immunoglobulin peptides in colon cancer. J. Proteome Res. 2015, 14 (9), 3555−3567. (13) Sheynkman, G. M.; Shortreed, M. R.; Frey, B. L.; Scalf, M.; Smith, L. M. Large-scale mass spectrometric detection of variant peptides resulting from nonsynonymous nucleotide differences. J. Proteome Res. 2014, 13 (1), 228−240. (14) Karpova, M. A.; Karpov, D. S.; Ivanov, M. V.; Pyatnitskiy, M. A.; Chernobrovkin, A. L.; Lobas, A. A.; Lisitsa, A. V.; Archakov, A. I.; Gorshkov, M. V.; Moshkovskii, S. A. Exome-driven characterization of the cancer cell lines at the proteome level: the NCI-60 case study. J. Proteome Res. 2014, 13 (12), 5551−5560. (15) Cox, J.; Hein, M. Y.; Luber, C. A.; Paron, I.; Nagaraj, N.; Mann, M. Accurate proteome-wide label-free quantification by delayed normalization and maximal peptide ratio extraction, termed MaxLFQ. Mol. Cell. Proteomics 2014, 13 (9), 2513−2526. (16) Craig, R.; Beavis, R. C. TANDEM: matching proteins with tandem mass spectra. Bioinformatics 2004, 20 (9), 1466−1467. (17) Kim, S.; Pevzner, P. A. MS-GF+ makes progress towards a universal database search tool for proteomics. Nat. Commun. 2014, 5, 5277. (18) Ivanov, M. V.; Levitsky, L. I.; Tarasova, I. A.; Gorshkov, M. V. Pepxmltka format converter for peptide identification results obtained from tandem mass spectrometry data using X!Tandem search engine. J. Anal. Chem. 2015, 70 (13), 1598−1599. (19) Jeong, K.; Kim, S.; Bandeira, N. False discovery rates in spectral identification. BMC Bioinf. 2012, 13 (Suppl 16), S2. (20) Goloborodko, A. A.; Levitsky, L. I.; Ivanov, M. V.; Gorshkov, M. V. Pyteomics–a Python framework for exploratory data analysis and rapid software prototyping in proteomics. J. Am. Soc. Mass Spectrom. 2013, 24 (2), 301−304. (21) Elias, J. E.; Gygi, S. P. Target-decoy search strategy for increased confidence in large-scale protein identifications by mass spectrometry. Nat. Methods 2007, 4 (3), 207−214. (22) Ivanov, M. V.; Levitsky, L. I.; Lobas, A. A.; Panic, T.; Laskay, Ü . A.; Mitulovic, G.; Schmid, R.; Pridatchenko, M. L.; Tsybin, Y. O.; Gorshkov, M. V. Empirical multidimensional space for scoring peptide spectrum matches in shotgun proteomics. J. Proteome Res. 2014, 13 (4), 1911−1920. (23) Käll, L.; Canterbury, J. D.; Weston, J.; Noble, W. S.; MacCoss, M. J. Semi-supervised learning for peptide identification from shotgun proteomics datasets. Nat. Methods 2007, 4 (11), 923−925.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel/Fax: +7-499-1378257. ORCID

Mikhail V. Gorshkov: 0000-0001-9572-3452 Notes

The authors declare no competing financial interest. All tables with filtered PSMs and peptides and X!Tandem output pepXML files are available at http://pubdata. theorchromo.ru/variantfdrcontrol/.



ACKNOWLEDGMENTS The work was supported by Russian Foundation for Basic Research (project #16-54-21006). S.M. and D.K. also thank the Russian Federal Agency for Scientific Organizations for the infrastructure support.



REFERENCES

(1) Zhang, B.; Wang, J.; Wang, X.; Zhu, J.; Liu, Q.; Shi, Z.; Chambers, M. C.; Zimmerman, L. J.; Shaddox, K. F.; Kim, S.; et al. Proteogenomic characterization of human colon and rectal cancer. Nature 2014, 513 (7518), 382−387. (2) Moghaddas Gholami, A.; Hahne, H.; Wu, Z.; Auer, F. J.; Meng, C.; Wilhelm, M.; Kuster, B. Global proteome analysis of the NCI-60 cell line panel. Cell Rep. 2013, 4 (3), 609−620. (3) Jaffe, J. D.; Berg, H. C.; Church, G. M. Proteogenomic mapping as a complementary method to perform genome annotation. Proteomics 2004, 4 (1), 59−77. (4) Kandoth, C.; McLellan, M. D.; Vandin, F.; Ye, K.; Niu, B.; Lu, C.; Xie, M.; Zhang, Q.; McMichael, J. F.; Wyczalkowski, M. A.; et al. Mutational landscape and significance across 12 major cancer types. Nature 2013, 502 (7471), 333−339. (5) Abaan, O. D.; Polley, E. C.; Davis, S. R.; Zhu, Y. J.; Bilke, S.; Walker, R. L.; Pineda, M.; Gindin, Y.; Jiang, Y.; Reinhold, W. C.; et al. The exomes of the NCI-60 panel: a genomic resource for cancer biology and systems pharmacology. Cancer Res. 2013, 73 (14), 4372− 4382. (6) Huang, P.-J.; Lee, C.-C.; Tan, B. C.-M.; Yeh, Y.-M.; Julie Chu, L.; Chen, T.-W.; Chang, K.-P.; Lee, C.-Y.; Gan, R.-C.; Liu, H.; et al. CMPD: cancer mutant proteome database. Nucleic Acids Res. 2015, 43 (D1), D849−D855. (7) Woo, S.; Cha, S. W.; Na, S.; Guest, C.; Liu, T.; Smith, R. D.; Rodland, K. D.; Payne, S.; Bafna, V. Proteogenomic strategies for identification of aberrant cancer peptides using large-scale nextgeneration sequencing data. Proteomics 2014, 14 (23−24), 2719− 2730. 1943

DOI: 10.1021/acs.jproteome.6b01014 J. Proteome Res. 2017, 16, 1936−1943