PeptideMapper - Efficient Amino Acid Sequence Mapping
PeptideMapper is a lightweight and efficient application for the mapping of amino acid sequences onto protein databases in the FASTA format. It takes advantage of the latest advances in string pattern matching to create a self-contained index that can be rapidly queried.
We have implemented queries for peptide sequences as well as sequence tags (alternating series of amino acids and mass gaps as produced by de novo sequencing algorithms).
PeptideMapper can be used standalone in command line or integrated as a dependency in other applications.
When using PeptideMapper, please cite Kopczynski et al., Bioinformatics, 2017.
For command line options and format specifications, please refer to PeptideMapperCLI. For use as a dependency, please refer to Usage As Dependency.
Methods
Database Index Creation
The first step is to prepare the protein sequences by creating one string with the pattern: /protein1/protein2/…/proteinn/ where proteini is the ith protein sequence in the given database and ‘/’ a separation character preventing the mapping of peptides on two consecutive proteins. Accordingly, we compute the suffix array (SA) [1] of the complete proteome using divsufsort, chosen for its performance. Subsequently, the Burrows and Wheeler transform (BWT) [2] is computed and stored in a wavelet tree (WT) data structure [3] for fast rank queries.
Only the 22 amino acid characters, the four characters for amino acid combinations and ‘/’ complete the alphabet Ε for the WT. To save memory, the WT is constructed according to a corresponding entropy-coded tree [4] built from the number of occurrences of every character. Additionally, the accession numbers extracted from the header are stored for every protein.
We implemented an efficient function querying all character occurrences within a substring in O(σ) where σ_=|Ε|. The SA is necessary to retrieve the starting positions of the mapped peptides within a protein. To save memory, we sample the SA and can retrieve the missing entries by using the last-to-front (LF)-mapping function [5] on the BWT.
For both peptide and sequence tag mapping, we utilize the backward search [5] which finds all occurrences of a pattern P in a text T in O(p) where p = |P|.
Peptide Mapping
To manage ambiguous amino acids and amino acids with indistinguishable mass, we first split peptides into their constituent amino acids and create a set of possible amino acids at each position. E.g., ATDWIRK => {A,X}{T,X}{D,B,X}{W,X}{I,L,J,X}{R,X}{K,X}. Theoretically, for this short peptide, we would need to map all possible 384 combinations, but using a dynamic programming approach we omit mapping of combinations containing prefixes we previously discarded.
The result is a set of ranges within the suffix array indicating matches that can be back calculated to the actual positions in the corresponding proteins. Having the position, we retrieve the accession number that completes all necessary data for further processing.
Sequence Tag Mapping
Sequence Tags are defined as alternating series of masses and sequences of amino acids. They can be inferred from spectra using sequence tag algorithms or de novo sequencing tools. A simple example is a tag of the following pattern m1-s-m2, where mi is a mass and s a partial peptide sequence. Having mapped s as described previously, we try to map the mass m1. Therefore, we extend s with all possible preceding sequences in the database creating a set of sequences s’.
For every sequence in s’, we store a nominal mass based on the individual amino acid masses and fixed modifications provided. We compute alternative possible masses based on the variable modifications. A sequence in s’ is kept when reaching m1 making a new set of sequences s”, and discarded if the nominal mass exceeds m1 plus a given tolerance mass. Subsequently, we extend every sequence in s” with its possible following amino acids as done previously until reaching m2. Since this extension is conducted in the opposite direction, we also compute an index for the reversed protein sequences. Having extended s to all sequences matching the tag, all proteins it appears in as well as the corresponding starting positions are returned.
Performance
For the following benchmark tests, we are using an Intel(R) Xeon(R) 2.80 GHz quad core desktop computer with 16 GB RAM. Note that only one core is used for the performance tests, since the task of mapping independent peptides against a sequence database can be heavily parallelized, the computation time can almost be divided by the number of cores used. We compared both compomics-utilities index methods, the ProteinTree and the new adapted FM-Index. Additionally, we compared the sequence tag mapping with TagRecon [6].
Protein Sequences Databases and Benchmark Datasets
Yeast (May 2016):
- size: 6,721 Proteins, 3,025,143 residues, 0 X residues
- fixed modifications: none
- variable modifications: none
- fragment tolerance: 0.02 Da
Mouse (May 2016):
- size: 16,813 Proteins, 9,474,758 residues, 79 X residues
- fixed modifications: Carbamidomethylation of C
- variable modifications: none
- fragment tolerance: 0.02 Da
Human (July 2015):
- size: 20,207 Proteins, 11,326,153 residues, 670 X residues
- fixed modifications: Carbamidomethylation of C, Oxidation of M
- variable modifications: none
- fragment tolerance: 0.02 Da
Proteogenomics (March 2015) [7], [8]:
- size: 83,721 Proteins, 13,851,427 residues, 0 X residues
- fixed modifications: Carbamidomethylation of C, Oxidation of M, Acetylation of K
- variable modifications: Acetylation of peptide N-term, Pyrolidone from Q, Pyrolidone from E
- fragment tolerance: 0.5 Da
Metaproteomics (January 2013) [9]:
- size 55,152 Proteins, 100,955,085 residues, 2,561,698 X residues
- fixed modifications: Carbamidomethylation of C
- variable modifications: Oxidation of M, Pyrolidone from Q, Acetylation of peptide N-term
- fragment tolerance: 0.2 Da
All UniProt proteomes concatenated (July 2016):
- size: 551,705 Proteins, 197,114,987 residues, 8,027 X residues
- fixed modifications: Carbamidomethylation of C, Oxidation of M
- variable modifications: Acetylation of K
- fragment tolerance: 0.02 Da
| For every database listed above, benchmark data sets D = {D1, D2, D3, D4} of different sizes ( | D1 | = 1,000, | D2 | = 10,000, | D3 | = 100,000, | D4 | = 1,000,000) were created by extracting peptides and tags at random positions in the database. All benchmark data (proteomes, data sets, parameter files) are available here. |
Results - Index Creation
| Database | time [s] | size [MB] | ||
|---|---|---|---|---|
| ProteinTree | FM-Index | ProteinTree | FM-Index | |
| Yeast | 65.73 | 2.55 | 184.3 | 7.02 |
| Mouse | 212.9 | 6.57 | 601.5 | 22.02 |
| Human | 252.0 | 7.49 | 537.9 | 26.33 |
| Proteogenomics | 526.2 | 8.81 | 1211 | 32.38 |
| Metaproteomics | 2956 | 58.00 | 5423 | 238.0 |
| All Proteomes | 7924 | 122.0 | 8653 | 458.0 |
Results - Peptide Sequences Mapping Time
| Database | D1 [s] | D2 [s] | D3 [s] | D4 [s] | ||||
|---|---|---|---|---|---|---|---|---|
| ProteinTree | FM-Index | ProteinTree | FM-Index | ProteinTree | FM-Index | ProteinTree | FM-Index | |
| Yeast | 6.067 | 0.041 | 16.96 | 0.206 | 22.45 | 1.385 | 72.01 | 12.01 |
| Mouse | 10.17 | 0.057 | 22.48 | 0.313 | 35.93 | 2.291 | 90.23 | 19.49 |
| Human | 8.260 | 0.056 | 25.85 | 0.403 | 123.1 | 2.297 | o.o.m. | 21.35 |
| Proteogenomics | 8.910 | 0.058 | 39.26 | 0.262 | 259.2 | 2.094 | o.o.m. | 17.71 |
| Metaproteomics | 70.72 | 0.224 | 1183 | 1.547 | 9015 | 13.49 | o.o.m. | 134.5 |
| All Proteomes | o.o.m. | 0.119 | o.o.m. | 0.643 | o.o.m. | 6.270 | o.o.m. | 62.20 |
Results - Sequence Tags Mapping Time
| Database | D1 [s] | D2 [s] | D3 [s] | D4 [s] | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| TagRecon | ProteinTree | FM-Index | TagRecon | ProteinTree | FM-Index | TagRecon | ProteinTree | FM-Index | TagRecon | ProteinTree | FM-Index | |
| Yeast | 238.0 | 4.670 | 0.096 | 279.0 | 20.64 | 0.590 | 281.0 | 38.12 | 3.863 | 288.0 | 205.9 | 33.59 |
| Mouse | 915.0 | 7.140 | 0.167 | 874.0 | 29.48 | 1.196 | 1002 | 107.3 | 9.910 | 1226 | 871.5 | 97.10 |
| Human | 1063 | 8.030 | 0.178 | 1151 | 31.76 | 1.352 | 1328 | 146.0 | 11.67 | 1608 | o.o.m. | 114.4 |
| Proteogenomics | PTMs not supported | 9.850 | 0.297 | PTMs not supported | 40.11 | 1.546 | PTMs not supported | 234.7 | 12.16 | PTMs not supported | o.o.m. | 119.1 |
| Metaproteomics | PTMs not supported | 63.93 | 6.281 | PTMs not supported | o.o.m. | 61.37 | PTMs not supported | o.o.m. | 609.3 | PTMs not supported | o.o.m. | 6205 |
| All Proteomes | 21981 | o.o.m. | 1.417 | 37791 | o.o.m. | 12.85 | 37353 | o.o.m. | 112.8 | 55171 | o.o.m. | 1145 |
o.o.m. = out of memory
Quick Start
To run PeptideMapper, download the packed zip archive of the latest compomics-utilities library, version 5.0.39 or newer.
wget http://genesis.ugent.be/maven2/com/compomics/utilities/5.0.39/utilities-5.0.39.zip
unzip utilities-5.0.39.zip
cd utilities-5.0.39
java -cp utilities-5.0.39.jar com.compomics.cli.peptide_mapper.PeptideMapperCLI -p exampleFiles/PeptideMapping/yeast.fasta exampleFiles/PeptideMapping/yeast-pep-1k.csv results.csv
A detailed description of the command line instructions is available here.
Troubleshooting
Should you encounter any issue with the usage of PeptideMapper, please create a new entry in the issue tracker.
References
[1] Manber and Myers, 1st Annual ACM-SIAM Symposium on Discrete Algorithms, 1990
[2] Burrows and Wheeler, Technical report, Digital Equipment
Corporation, 1994
[3] Grossi et al., Proceedings of the Fourteenth Annual
ACM-SIAM Symposium on Discrete Algorithms, 2003
[4] Huffman, Proceedings of the Institute of Radio Engineers, 1952
[5] Ferragina and Manzini, Proceedings of the 41st Annual Symposium on Foundations of Computer Science, 2000
[6] Dasari et al., Journal of Proteome Research, 2010
[7] Menschaert et al., Molecular & Cellular Proteomics, 2013
[8] Crappé et al., Nucleic Acids Research, 2014
[9] Tanca et al., PLoS ONE, 2013