Genome-based peptide fingerprint scanning

Computational Methods

In silico digestion of an entire genome sequence generates masses for all peptides that might be produced by its in vivo translation and subsequent proteolytic digestion. Translation and digestion are performed in all forward and reverse frames. The process proceeds from 5′ to 3′ on each sequence, keeping an in-process list of not yet terminated fragments and a final master list of terminated fragments. As each in-frame codon is encountered the mass of its amino acid translation is added to each fragment of the in-process list. A new fragment is added to the incomplete list after each cleavage or stop codon and at each start codon. Terminated fragments are transferred to the master list when cleavage sites and stop codons are encountered. Any fragment falling below a length threshold of three codons is discarded.

The enzyme trypsin cleaves after lysine and arginine, encoded by a total of eight codons. In the standard genetic code, the start codon is ATG (also coding internal methionine) and stop codons are TAA, TAG, and TGA. For mtDNA, start and stop codons are ATG, ATA and TAA, TAG, respectively. The system supports the use of alternative genetic codes by using separate translation dictionaries, selected by annotation in the FASTA header of a sequence, and defaulting to standard nuclear encoding.

Duplicate fragments (produced, e.g., when a start codon follows a cleavage codon) are prevented by checking a queue containing fragments recently transferred for any fragment with the same start codon and mass before transferring the new fragment to the final digest list. Incomplete proteolytic digestion of a protein sample results in peptides containing internal cleavage sites. The presence of such missed cleavages requires that the in silico digest peptides also contain them. The program uses the variable b to specify the maximal number of missed cleavage sites (breaks) allowed internally to a fragment. For efficiency, b is generally kept low, i.e., b ≤ 2.

In silico digestion of the nuclear plus mitochondrial genome of S. cerevisiae with b = 2 generates 8.9 million fragments and requires ≈200 MB of memory. The program stores the mass, start point, length, number of breaks, and reading frame of each fragment. The computational complexity of the digestion process is roughly linear, or O(n), in correspondence to genome size n. The matching process is O(nm), with m representing the size of the mass list. Further algorithmic details are provided in Supporting Text, which is published as supporting information on the PNAS web site, www.pnas.org.

Mass data from either MALDI-MS or ESI-MS is matched to the in silico digestion to within a tolerance level calculated as a percentage of the experimental mass (Δ%). Each matching fragment is mapped back to its position on the chromosome (Fig. 1). Two scans across each chromosome evaluate match criteria for windows of size w to identify those with a high score (scoring discussed below). For all experiments herein w = 500. The first scan scores windows in 100-nt increments, providing an internal histogram of score distributions used by the program to establish a cutoff level such that only the top n scoring clusters are examined on the subsequent scan. Currently, n is fixed at a value of 10. Regions with windows scoring above the threshold are examined in further detail by an extension scan used to determine the full extent of the hit-cluster region. The extension scan starts with each high-scoring window and proceeds backward in 50-nt steps, considering each window of size w until the score falls below a defined cutoff. This marks the start of the full region. A scan forward from the original window in 50-nt steps again considers windows of size w until the same cutoff is reached, marking the end of the extended region. The cutoff currently used is half of the score of the 10th highest-scoring window found on the initial scan.

Several scoring methods have been investigated. Simple measures include the total number of hits (matches) per window and the percentage of the DNA sequence contained within matching fragments. The former tends to favor regions containing multiple repeats of a peptide matching one of the input masses and fails to account for the lower frequency of peptides with missed cleavages. Sequence coverage scores are inflated by high-mass fragments. A scoring function was developed to address these issues. It considers the following aspects of each window of size w: (i) the number of hits in a single frame; (ii) the number of possible hits in the same frame (i.e., the total digested fragment count); (iii) the number of missed cleavages in each fragment; (iv) the number of in-frame stop codons encountered in the window before the current fragment; (v) duplicate mass matches; and (vi) abutment of fragments.

Multiplicative combination of these attributes is important for scoring features such as the number of preceding stop codons, missed cleavages, and duplicate fragments matched. For example, whereas a region with one in-frame stop may still be considered because it could be caused by a sequence error or stop codon read-through (13, 20), multiple in-frame stops are increasingly unlikely. If the probability of stop-codon present in-frame is p, then the probability of s occurrences is p^s. The scoring equation does not attempt to model actual probabilities. The data required to ascertain realistic probabilities or frequencies of these occurrences would be very difficult to obtain. Instead, the scoring function is intended to maximally discriminate randomly formed clusters from real hits, while allowing for occurrences such as sequence errors or missed cleavages.

Assignment of penalty values to each of the listed factors allows their multiplicative combination. The values are: c_b, penalty for missed cleavages, default 0.6; c_s, penalty for preceding in-frame stops, default 0.4; c_d, penalty for preceding duplicate-mass matches, default 0.6; and c_ā, penalty for N terminus not abutting a preceding fragment, default 0.9.

For a window containing t hypothetical fragments, h of which match experimentally measured peptide masses, we calculate a window score s by summing the penalty products for each fragment j = 1 … h: 1 The functions b(f_j), s(f_j), d(f_j), and a(f_j) return counts for the number of breaks in a fragment (e.g., Lys and Arg codons), the number of stops preceding a fragment, the number of other preceding matches for this mass in the window, and whether or not (1 or 0) the amino terminus of the fragment abuts a preceding fragment, respectively. The term t normalizes the results for the total number of digest fragments in the window, but used by itself can skew results toward windows with a small number of possible fragments; h counterbalances this by giving weight to the number of hits in the window, and is reduced by d, the number of duplicate matches in the window. Scores are multiplied by a scaling factor of 100 to simplify the histogram analysis.

These scores are used to differentiate statistically significant match regions from the backdrop of random hits. Methods to assign probability estimators to fingerprint matches have been described by several groups. One approach is to derive a probabilistic function describing the likelihoods of alternative identifications, as was done in profound (5), or to assign a probability that a given protein match is by chance as was done in mascot (6). Another is to use randomized data to establish a baseline for assigning significance to matches with real data (21).

We use a method similar to the latter, by which we calculate the significance of a match region as a function of its window score and the total number of masses in the experimental spectrum. Each value in a randomly selected subset of a large set of experimental peptide masses is perturbed by a random modification representing the addition or subtraction of H, C, O, and N atoms. A large number of such mass lists is searched against the genome fragment database, establishing a histogram to represent the range of the null hypothesis (i.e., that any given result is caused by chance). The histogram is used to define P_s, the probability in a single genome scan that a randomly chosen set of masses, of the same size, would achieve a score equal or above the score considered (derivation in Supporting Text). The scoring methods used by GFS and mascot cannot be directly equated because the GFS significance score describes the probability of a false positive in a complete genome scan, whereas mascot produces a score that describes the probability of false-positive for each protein considered.

The system consists of a client-server pair with a UNIX command-line interface. The server performs the in silico digestion and keeps the resulting database in memory, obviating the recomputation of the genome digest over multiple MS analyses. Each client receives a peptide mass list, connects to the server for processing, and outputs the results. The client generates a formatted HTML file displaying the clusters for the 10 highest-scoring windows, with matched fragments highlighted in different colors according to the reading frame in which they are found to facilitate visual identification of transframe events. The entire region's score and the highest score for any contained window of size w are reported. The significance is calculated as a function of the latter, maximal fixed window-size score.

Server parameters include the number of missed cleavages to be calculated, window size, a directory containing all FASTA-formatted genome sequences, and whether to use average or mono-isotopic masses. Client parameters include the file containing a peptide mass list, the host name of the server, the tcp port number, Δ% mass tolerance, and parameters related to random trials. The programs currently run on the UNIX command line of MacOS X, with potential deployment on Linux. We plan to make a graphical interface accessible to other researchers via the web. The current prototype code is available free of charge to nonprofit researchers.

Laboratory Methods

Results

Data from ESI-MS analysis of S. cerevisiae mitochondrial proteins and MALDI-MS analysis of E. coli proteins were used to test the GFS system. We rejected the use of synthetic data for reasons similar to those cited by Perkins et al. (6), e.g., because the real experimental factors that play into the observed data are not well enough understood to be modeled as required for the generation of realistic synthetic data.

On the other hand, the absence of authoritative identifications for the analyzed proteins requires that we evaluate performance by comparing our results with those of other algorithms. Our performance assessment is based on comparisons with both peptident (4) and mascot (6), well recognized tools used for the data analysis in our yeast proteome project. In most cases both of the programs were used for comparison; however, there are a few for which only mascot was used, because later in the S. cerevisiae project that became the default analysis tool. For mascot analyses in our proteome project, proteins with a match score >70 were considered significant identifications.

As opposed to the single protein samples typically produced by 2D gel separation, our yeast mitochondrial samples contain multiple proteins, a result of the method of HPLC separation used. The regular presence of multiple proteins imposes certain constraints on the system. For example, it greatly complicates the consideration of mutually exclusive possibilities required to develop a Bayesian formula for probability assessment of identifications as was done in profound (5). An alternative is to base comparisons on the distribution of randomized data, as we have done here. The increased number of peptides in a multiprotein analyte further increases the background noise level, making matches more difficult to distinguish.

To establish the significance levels of scores we performed multiple repetitions of experiments by using randomized mass lists of varying lengths. For S. cerevisiae, 1,000 repeat trials were performed, with lists of length 20–100 in increments of five. An example histogram of scores produced from one such set of trials is shown in Fig. 3, which is published as supporting information on the PNAS web site. The process was repeated for E. coli with lists of 20–50 random masses, incremented by five. The program keeps a summary histogram plotting scores against the number of w = 500-nt windows achieving each score. A plot of the P_s< 0.001 (maximum scores) and P_s< 0.05 values versus mass-list size is in Fig. 2. The curves show a close to linear relationship between the number of masses input and confidence thresholds, with the P_s< 0.05 varying more smoothly because of the higher quantity of data available to establish it.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Confidence levels established by sets of random mass-list trials against the S. cervisiae whole-genome digest. ⧫, Top scores in 1,000 trials, equating to P_s < 0.001. ■, The scores above which P_s < 0.05. The dotted and solid lines are a linear regression for the two data sets, shown along with their equations and R values.

Twenty-two samples were chosen for the comparative performance analysis, 18 from yeast and four from E. coli. The yeast samples were selected at random from the larger pool of those available. The program was run with each mass list, and results were manually parsed, with the top two scoring clusters considered in each case. The score for each cluster region was compared with the null hypothesis to compute the P value by interpolating between the two nearest mass-list sizes (increments of five).

The genome position of each identified cluster was compared with ORF annotations for yeast or E. coli. An ORF encompassing and in-frame with the cluster region was noted as a match. For the yeast searches, an in-house program was used to match positions to annotated ORFs in a mid-2001 download from the Saccharomyces Genome Database (http://genome-www.stanford.edu/Saccharomyces/). E. coli searches were performed manually against the National Center for Biotechnology Information database by using the summary data at www.ncbi.nlm.nih.gov/cgi-bin/Entrez/altik?gi=115&db=Genome (late 2001). The parameters for all reported experiments were w = 500, Δ% = 0.05, and b = 3, with monoisotopic peptide masses used for the yeast/ESI experiments and average peptide masses for the E. coli/MALDI experiments. Fig. 4, which is published as supporting information on the PNAS web site, is an illustration of typical output from the program.

Seventeen of the 22 samples had a top-scoring cluster region with significance P_s< 0.05. Of the five not within this threshold, four had top-scoring clusters falling in ORFs identified in the database as mitochondrially localized proteins; the remaining region did not correspond to an obvious ORF. The top-scoring GFS-identified ORF was also top scoring for one of the other programs in 16 cases, including one case where mascot and GFS agreed that no significant matches were present. When counting agreement between either first- or second-place significant protein matches, the GFS and mascot/peptident corroboration increases to 19 samples or 86%. Table 1 shows six representative results, including two disagreements and one where the second-place identification was corroborated. The disagreements provide insights into differences between the algorithms. Sample 2378 F9 with GFS score 121 (P_s ≈ 0.004) mapped to an ORF for the mitochondrial precursor of alcohol dehydrogenase (ADH3), whereas the other programs identified the sample protein as phosphoglycerate kinase (PGK). Given the high significance and mitochondrial localization of ADH3, it is likely that GFS correctly identified a different component of the sample than the other programs. GFS also corroborated the PGK identification, it being the second top-scoring region found.

Another analyte was identified by GFS as PET127, a protein annotated in GenBank as a component of the mitochondrial translation system. For this sample neither mascot nor peptident found anything significant. It appears from the GFS program output that only a small, ≈10-kDa portion of the protein was present during tryptic digestion, which would account for the poor significance score. It would also explain the difficulty mascot had with this, because this is a small piece of a large 93-kDa protein predicted in the ORF database. In another case, GFS identified POM152 (P_s ≈ 0.03), whereas mascot found nothing significant and peptident had a weak match for ABC1. As with the previous case, the parent protein predicted by the database (Saccharomyces Genome Database) is large: 151 kDa. This finding highlights an advantage of position-based peptide matching. It is likely that a portion of analyzed proteins has undergone in vivo proteolysis, causing incomplete peptide coverage. Because the peptide coverage, when averaged over a large protein, will be low, such cases can confound searches that rely on ORF or protein annotation. With the genome-based positional scanning of GFS, parent protein size does not directly affect its performance unless a protein is much smaller in size than the window used for analysis.

Table 2 illustrates the detection of multiple proteins within a sample containing a large number of peptides (88 total). To calculate the significance for each match region we remove from the experimental mass list all masses matched to a higher-scoring region and not contained in the current region and calculate the significance corresponding to the score and this new number of masses. The rationale for this is its equivalency to removing the masses previously matched and rerunning the scan. GFS is unique among MS search algorithms for its ability to establish significance levels based only on the number of masses input. The use of a fixed window size for scanning simplifies the detection of multiple proteins and assessment of their significance in a single-pass analysis.

Identifications for all four E. coli samples matched the standard results, and all of them had GFS scores with significance P_s< 0.001. The stronger E. coli results are likely caused by the samples each having only one protein and by higher mass accuracy from both the instrument and data analysis methods. The ESI process generates ions with multiple charge states, whereas MALDI typically produces singly charged ions. For the ESI data, we lacked appropriate software to perform deconvolution of multiply charged spectra into straight mass spectra; our manual deconvolution may have been less accurate. Given the lower quality of these data and the relaxed match criteria used (500 ppm), the performance of GFS on the ESI/yeast data is promising.

Discussion

The empirically determined working parameter sets sufficed for practical operation on these distinct data sets but should be optimized by further experimentation to match the properties of the experimental equipment used. For example, mass tolerance should be decreased for higher mass accuracy data, effectively reducing the number of random hits. Optimization of window size is potentially complex, depending on protein size, gene composition, and the local configuration of matched fragments. If spliced genes are analyzed, the window size may have to be much larger, with a second, smaller window scan to identify putative exons.

We performed an experiment to investigate the extent of the random backdrop for large mammalian genomes. We repeated 1,000 queries with randomly generated 41- and 61-length mass lists against the in silico digest of human chromosome XIV, comprising 1/35th of the human genome. The background from this experiment should be roughly equivalent to scanning the entire genome 28 times. The maximum score for 41-mass lists was 107, and for 61-mass lists was 134, indicating that a real match scoring above those thresholds would have a significance value of less than ≈1/28 = 0.04. Our yeast sample (2378 E2) that had a 41-mass input list produced a score of 114, which for yeast is a P_s value <0.001, and for a genome of human size and composition would still have significance <0.04. A 61-mass sample that had a P_s≈ 0.004 for yeast would have a P_s≈ 0.13 if matched against a human-sized genome. These data indicate that scanning a much larger genome is statistically feasible. Improved MS data, allowing reduction of the tolerance from the present 500 ppm to ≈50 ppm, should provide a proportional 10-fold reduction in the random backdrop, improving the likelihood of success. However, the identification of proteins consisting of multiple exons remains a considerable challenge that has not yet been addressed.

Although our implementation is a proof-of-principle prototype, we have obtained interesting results contributing useful information to our research of yeast mitochondrial proteins. In application to a less-annotated genome than yeast, this method could contribute to the identification of proteins for which the database representation is not complete. Potential future enhancements include the automatic assignment of probability values from the regression of random trials, integration with ORF databases for automatic output of the ORF covering any high scoring cluster, and a web-based interface. We are optimizing the code for high-throughput analysis on modern vector processors and plan to deploy the system to provide concurrent search of multiple mass lists and several genomes.