Skip to content

PeptideMapper

Jump to bottom
dominik-kopczynski edited this page Jan 25, 2017 · 129 revisions

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.

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 index methods the prior ProteinTree with the new adapted FM-Index.

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) [6]:

  • 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) [7], [8]:

  • 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.35 5423 238.7
All Proteomes 7924 122.08 8653 458.1

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.018
Mouse 10.17 0.057 22.48 0.313 35.93 2.291 90.23 19.49
Human 8.26 0.056 25.85 0.403 123.1 2.297 o.o.m. 21.35
Proteogenomics 8.91 0.058 39.26 0.262 259.2 2.094 o.o.m. 17.71
Metaproteomics 70.72 0.224 1183.8 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]
ProteinTree TagRecon FM-Index ProteinTree TagRecon FM-Index ProteinTree TagRecon FM-Index ProteinTree TagRecon FM-Index
Yeast 4.67 238 0.096 20.64 279 0.59 38.12 281 3.863 205.9 288 33.596
Mouse 7.14 915 0.167 29.48 874 1.196 107.3 1002 9.910 871.5 1226 97.10
Human 8.03 1063 0.178 31.76 1151 1.352 146.0 1328 11.67 o.o.m. 1608 114.4
Proteogenomics 9.85 PTMs not supported 0.297 40.11 PTMs not supported 1.546 234.7 PTMs not supported 12.16 o.o.m. PTMs not supported 119.1
Metaproteomics 63.93 PTMs not supported 6.281 o.o.m. PTMs not supported 61.37 o.o.m. PTMs not supported 609.3 o.o.m. PTMs not supported 6205
All Proteomes o.o.m. 21981 1.417 o.o.m. 37791 12.85 o.o.m. 37353 112.8 o.o.m. 55171 1145.0
o.o.m. = out of memory

Quick run

To run PeptideMapper, download the packed zip archive of the latest compomics-utilities library, at least 4.9.0.

wget http://genesis.ugent.be/maven2/com/compomics/utilities/4.9.0/utilities-4.9.0.jar
unzip utilities-4.9.0.jar
cd utilities-4.9.0
java -cp utilities-4.9.0.jar com.compomics.util.experiment.identification.protein_inference.executable.PeptideMapping -p exampleFiles/PeptideMapping/yeast.fasta exampleFiles/PeptideMapping/yeast-pep-1k.csv results.csv

A detailed description of the comand 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] Tanca et al, PLoS ONE, 2013
[7] Menschaert et al, Molecular & Cellular Proteomics, 2013
[8] Crappé et al, Nucleic Acids Research, 2014

Clone this wiki locally