PepQuery enables fast, accurate, and convenient proteomic validation of novel genomic alterations

PepQuery workflow

PepQuery takes as input a novel peptide, protein, or DNA sequence, or novel genomic features in the VCF, BED, or GTF file format, and the workflow includes five major steps: (1) target peptide sequence preparation and initial filtering; (2) candidate spectra retrieval and PSM scoring; (3) competitive filtering based on reference sequences; (4) statistical evaluation; and (5) competitive filtering based on unrestricted modification searching (Fig. 2; Methods). For PSM scoring, we implemented the Hyperscore used in X!Tandem (Craig and Beavis 2004) and the multivariate hypergeometric distribution (MVH) score used in MyriMatch (Tabb et al. 2007). For statistical evaluation of the PSM scores, permutation P-values are calculated based on randomly shuffled peptide sequences (Methods). The spectra with statistically significant matches to the query peptide that cannot be better explained by unmodified or modified reference sequences are reported. Annotated PSMs can be visualized in the tool for manual evaluation.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

PepQuery workflow. (A) The PepQuery workflow involves five major steps: (1) target peptide sequence preparation and initial filtering; (2) candidate spectra retrieval and PSM scoring; (3) competitive filtering based on reference sequences; (4) statistical evaluation; and (5) competitive filtering based on unrestricted post-translational modification searching. The red text illustrates a real example in which a variant peptide LVVVGADGVGK is used to query the CPTAC colorectal cancer (CRC) data set with 95 samples and 12,941,421 spectra, using RefSeq as the reference protein database. Whereas existing methods require pairwise analysis between all 12,941,421 spectra and all RefSeq-derived peptide sequences plus the variant peptide sequence, PepQuery focuses only on the spectra that are relevant to the novel peptide, which reduces computational time and also allows more comprehensive analysis of these spectra. Illustration of peptide (B) and modification (C) indexing methods.

With the existing customized database approach, identifying a novel peptide sequence as shown in Figure 2 in the Clinical Proteomic Tumor Analysis Consortium (CPTAC) colorectal cancer data set (Zhang et al. 2014) will require constructing a customized database that includes both reference protein sequences (e.g., from RefSeq) and the novel sequence, and then evaluating all possible spectrum-peptide pairs between over 12 million MS/MS spectra and over one million unique peptide sequences. In contrast, PepQuery starts with the novel peptide sequence of interest, evaluates only peptide-spectrum pairs involving the novel peptide, and thus reduces the computational time by several orders of magnitude (Fig. 2). Moreover, the peptide-centric approach quickly narrows down to a small number of candidate spectra, which makes it possible to perform unrestricted modification searching in the last step to identify alternative interpretations of the candidate spectra in real-time applications.

Validating completely novel protein sequences

To evaluate the performance of PepQuery for validating completely novel protein sequences, we used a spiked-in MS/MS proteomic data set from the Proteome Informatics Research Group (iPRG) 2015 study (Choi et al. 2017). This data set was generated from yeast cells with six spiked-in nonyeast proteins. Here, we used the yeast proteome as the background, and the peptides derived from the six spiked-in nonyeast proteins were treated as the ground truth of “novel peptides.” The goal was to investigate whether PepQuery could accurately identify the “novel peptides” using individual “novel peptides,” the MS/MS data, and the yeast reference proteome as input. Because some of the peptides derived from the spike-in proteins may not be detectable in the MS/MS experiment, we first performed a standard database search using IPeak (Wen et al. 2015) to determine which peptides from the spiked-in proteins were detectable in the MS/MS experiment with 1% false discovery rate (FDR) at both PSM and protein levels. IPeak integrates search results from three search engines, MS-GF+, MyriMatch, and X!Tandem. Although this does not necessarily produce the exact set of detectable peptides from the spike-in proteins, using multiple search engines has been shown to achieve higher sensitivity and specificity compared to individual search engines (Wen et al. 2015). The IPeak analysis identified a total of 93 peptides from the six spiked-in proteins (Supplemental Table S1). The 93 peptides were used as gold-standard positives for novel protein sequences (Supplemental Fig. S1). Meanwhile, we randomly selected 10,000 tryptic peptides from Escherichia coli as gold-standard negatives. Searching the gold-standard peptide sequences against the iPRG data set using PepQuery, we obtained true positive rates (TPRs, a.k.a., sensitivity) of 94.62% and 97.85% (Fig. 3A; Supplemental Tables S2, S3) and false positive rates (FPRs, a.k.a., 1-specificity) of 0.05% and 0.05% (Fig. 3B; Supplemental Tables S4, S5) based on the Hyperscore and MVH score, respectively. The precisions of using Hyperscore and MVH after unrestricted filtering were 94.62% and 94.79%, respectively. These results demonstrate high sensitivity and specificity of PepQuery.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Evaluation of the performance of PepQuery for validating completely novel protein sequences. (A) Sensitivity and (B) specificity evaluation for novel peptide sequence validation. Mod. filtering means unrestricted modification searching-based filtering. (C) A spectrum matched to a gold-standard negative peptide without unrestricted modification searching-based filtering. Blue and red colors indicate matched peaks, whereas gray indicates unmatched peaks. (D) The same spectrum in C can be matched to a reference peptide with an ammonium salt modification on aspartic acid with very few unmatched peaks and a high score.

Notably, filtering based on unrestricted modification searching reduced the false positives by 89% (from 47 peptides to five peptides) and 95% (from 93 peptides to five peptides) for the two scoring algorithms, respectively (Fig. 3B), whereas the true positives were only reduced by 1% (Fig. 3A). To illustrate the merit of this approach in reducing false positives, a spectrum incorrectly matched to an E. coli peptide (false positive) without unrestricted modification searching-based filtering (MVH score = 19.7, P = 0.0028) (Fig. 3C) was matched by PepQuery to a yeast peptide with an ammonium salt modification on aspartic acid (D) (MVH score = 70.0) (Fig. 3D). When including ammonium salt modification on aspartic acid as a variable modification in a standard database searching of the iPRG 2015 data set using MyriMatch, we found that 306 (1.33%) out of the 22,678 identified peptides had the same ammonium salt modification, suggesting its origin from artifacts of sample handling. All spectra incorrectly matched to an E. coli peptide without unrestricted modification searching-based filtering and the alterative matches to a modified yeast peptide were annotated and presented in Supplemental Figure S2. These data suggest that competitive filtering with unrestricted modification searching provides a powerful means to reduce false positives in proteogenomic studies.

Validating novel splice junctions

To evaluate the ability of PepQuery to validate novel splice junctions, we applied PepQuery to a data set from a published study on novel splice-junction peptide identification (Sheynkman et al. 2013). The MS/MS data were generated from the Jurkat cells, and the study reported a total of 55 unique splice-junction peptides (Supplemental Table S6). The 55 peptides were searched against the MS/MS data using PepQuery and 36 (65%) and 39 (71%) were confirmed based on the Hyperscore and MVH score, respectively (Fig. 4A; Supplemental Tables S7, S8). Among the 16 peptides not confirmed by MVH (i.e., false negatives), nine were removed by competitive filtering based on reference sequences and five were removed by competitive filtering based on unrestricted modification searching. Thus, 87% of the MVH-based “false negatives” were likely to be false positive discoveries in the original study (Fig. 4C) and so were 84% of the Hyperscore-based “false negatives” (Fig. 4B). These data demonstrate that PepQuery effectively removes false identifications reported in published proteogenomic studies.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Evaluation of the ability of PepQuery to validate novel splice junctions. (A) Novel splice-junction peptide identification result. (B,C) Classification of the novel junction peptides reported by the original study but not by PepQuery (orange sections in the large circles). I: Peptides removed by competitive filtering based on reference sequences. II: Peptides removed by competitive filtering based on unrestricted modification searching. III: Remaining false negatives.

Validating single amino acid variants

To evaluate the ability of PepQuery to validate single amino acid variants (SAAVs), we used simulation data derived from the iPRG 2015 study (Methods). For each of the 93 peptides from the six spiked-in nonyeast proteins, we generated two versions of variant peptides, each with a randomly introduced SAAV. One version served as a gold-standard negative. The other version was added to the reference database to replace the original reference peptide sequence. With regard to this new reference sequence, the original peptide became an SAAV peptide. Because the original peptide had matched spectra in the MS/MS data, it was used as a gold-standard positive sequence in our evaluation (Supplemental Fig. S3). On the basis of 100 sets of gold-standard positives, gold-standard negatives, and expanded reference databases, PepQuery achieved average TPRs of 93.53% and 96.27% (Fig. 5A) and average FPRs of 0.41% and 0.26% (Fig. 5B) using the Hyperscore and MVH score, respectively. Filtering based on unrestricted modification searching reduced the false positives by 75% and 86% for the two scoring algorithms, respectively (Fig. 5B), with only 1.56% and 1.96% reduction on the true positives, respectively (Fig. 5A). These results are comparable to those from the novel protein sequence study (Fig. 3A,B) and further highlight the sensitivity and specificity of PepQuery in SAAV validation and the benefits of competitive filtering with unrestricted modification searching.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Evaluation of the ability of PepQuery to validate single amino acid variants (SAAVs). (A) Sensitivity and (B) specificity evaluation for validating single amino acid variants based on 100 simulation studies.

Applying PepQuery to large proteomic data sets

An important emerging application of PepQuery is to identify proteomic evidence for genomic findings using large, publicly available proteomic data sets. To demonstrate the utility of PepQuery in this scenario, we used the colorectal cancer data set produced by CPTAC (Zhang et al. 2014). The 4084 SAAV peptide-sample pairs reported by the original paper, which involves 799 fully tryptic SAAV peptides and 79 samples, were used as gold-standard positives for this study, whereas all other 59,037 pairs between the 799 SAAV peptides and the 79 samples were used as gold-standard negatives (Methods; Supplemental Table S9). Using PepQuery, each of the 799 SAAV peptides was searched against all 10,753,601 MS/MS spectra from the 79 samples. Annotated spectra for PSMs identified based on Hyperscore and MVH score can be found in Supplemental Figures S4 and S5, respectively. The Hyperscore and MVH score-based analyses showed TPRs of 78.04% and 90.96% (Fig. 6A) and FDRs of 2.01% and 2.40% (Fig. 6B), respectively (Supplemental Tables S10, S11). The TPRs are lower than those observed in the simulation study, whereas the FPRs are higher.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

Evaluation of the performance of PepQuery for validating SAAVs in large data sets. (A) Sensitivity and (B) specificity evaluation for validating SAAVs in a large cancer data set. (C) A large proportion of the false positives in B are supported by RNA-seq evidence. (D,E) Two spectra supporting the KRAS G12D mutation in a sample, in which the mutation was supported by DNA but not RNA evidence.

We note that, unlike the simulation study, the “gold-standards” used in this analysis were derived from the original customized database searching results and are not perfect ground truth. In the original study, database searching used MyriMatch (MVH) but not Hyperscore, which explains the much lower TPR observed for the Hyperscore-based analysis. Furthermore, competitive filtering based on unrestricted modification searching reduced the true positives by 7.54% and 4.39% for the two algorithms, respectively (Fig. 6A). The original study considered only a few PTMs. Therefore, some of the gold-standard positives may be false because the spectra supporting the SAAV peptides corresponded better to modified reference sequences.

On the other hand, some of the false positives may be true. Among the SAAV peptide-sample pairs identified by PepQuery but not the original study (Fig. 6B), we were able to find RNA-seq evidence for 39% of the Hyperscore-based “false positives” (H1, 457 SAAV peptide-sample pairs) and 41% of the MVH score-based “false positives” (M1, 576 SAAV peptide-sample pairs) (Fig. 6C; Supplemental Tables S10, S11). These SAAV peptide-sample pairs were missed in the original report because RNA-seq data from these samples only had one or two reads covering corresponding SNVs, and a minimum of three-read depth was required for including an SNV in a sample-specific customized database (Zhang et al. 2014). We also found that unrestricted modification searching-based filtering increased the proportion of “false positives” that could be supported by RNA-seq evidence (Fig. 6C, H1 vs. H2 and M1 vs. M2).

RNA-seq also may fail to detect some SAAVs. For example, despite a lack of RNA-seq evidence, PepQuery identified two spectra supporting the KRAS G12D mutation in the sample TCGA-AA-A02O-01A (Fig. 6D,E). We manually checked the cBioPortal (Cerami et al. 2012), and the KRAS G12D mutation was indeed reported for this sample. Taken together, PepQuery identified hundreds of SAAV events that are supported by genomic data but are missing in the original customized database searching results, and unrestricted modification searching-based filtering helped remove false positives.

FDR estimation

To estimate FDR corresponding to the PepQuery P-value cutoff of 0.01 used in this study, we performed two additional analyses. The first analysis used a human-E. coli spike-in data set (317,502 MS/MS spectra) from a published study (Shen et al. 2017), in which E. coli protein lysate was spiked into a human cell digest. We searched the MS/MS data against the combined human-E. coli SWISS-PROT database using X!Tandem, MyriMatch, and MS-GF+ and identified 4540 unique E. coli peptides and 30,099 unique human peptides. The 4540 E. coli peptides were used as the target peptides. We further generated 4540 reverse peptides based on the target peptides. After removing reverse peptides that could be matched exactly to a human or E. coli protein sequence, we got 4512 decoy peptides. We searched the 9052 target and decoy peptide sequences against the human-E. coli spike-in data set using PepQuery and a P-value cutoff of 0.01. As shown in Table 1, the FDRs were 0.84% and 0.37% in the Hyperscore- and MVH score-based analyses, respectively.

In the second analysis, we further evaluated the impact of sample size on PepQuery FDR using the CPTAC colon cancer data set (Zhang et al. 2014) with a total of 79 samples. We searched the MS/MS data against the RefSeq human protein sequences database using X!Tandem, MyriMatch, and MS-GF+ and identified 77,069 unique human peptides. We performed FDR estimation for randomly selected sets of one sample, 19 samples, 39 samples, 59 samples, and 79 samples. For each sample set, we randomly selected 2000 target peptides and generated 2000 decoy peptides. We searched all 4000 peptides against the MS/MS data from the same sample set using PepQuery and a P-value cutoff of 0.01. The analysis was repeated five times for each sample set. We observed increased FDRs with increasing sample size number (Table 2). Nevertheless, the FDRs were still <3% when all 79 samples were analyzed. These results suggest that a 1% PepQuery P-value cutoff coupled with unrestricted modification searching-based filtering provides well-controlled FDR for small to relatively large sample sizes.