OpenMS  2.5.0
PeptideIndexing.h
Go to the documentation of this file.
1 // --------------------------------------------------------------------------
2 // OpenMS -- Open-Source Mass Spectrometry
3 // --------------------------------------------------------------------------
4 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
5 // ETH Zurich, and Freie Universitaet Berlin 2002-2020.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Chris Bielow $
32 // $Authors: Andreas Bertsch, Chris Bielow $
33 // --------------------------------------------------------------------------
34 
35 #pragma once
36 
37 
54 #include <OpenMS/SYSTEM/SysInfo.h>
55 
56 #include <atomic>
57 #include <algorithm>
58 #include <fstream>
59 
60 
61 namespace OpenMS
62 {
63 
124  class OPENMS_DLLAPI PeptideIndexing :
125  public DefaultParamHandler, public ProgressLogger
126  {
127 public:
128 
131  {
136  UNEXPECTED_RESULT
137  };
138 
140  PeptideIndexing();
141 
143  ~PeptideIndexing() override;
144 
145 
147  inline ExitCodes run(std::vector<FASTAFile::FASTAEntry>& proteins, std::vector<ProteinIdentification>& prot_ids, std::vector<PeptideIdentification>& pep_ids)
148  {
149  FASTAContainer<TFI_Vector> protein_container(proteins);
150  return run<TFI_Vector>(protein_container, prot_ids, pep_ids);
151  }
152 
188  template<typename T>
189  ExitCodes run(FASTAContainer<T>& proteins, std::vector<ProteinIdentification>& prot_ids, std::vector<PeptideIdentification>& pep_ids)
190  {
191  // no decoy string provided? try to deduce from data
192  if (decoy_string_.empty())
193  {
194  auto r = DecoyHelper::findDecoyString(proteins);
195  proteins.reset();
196  if (!r.success)
197  {
198  r.is_prefix = true;
199  r.name = "DECOY_";
200  OPENMS_LOG_WARN << "Unable to determine decoy string automatically (not enough decoys were detected)! Using default " << (r.is_prefix ? "prefix" : "suffix") << " decoy string '" << r.name << "'\n"
201  << "If you think that this is incorrect, please provide a decoy_string and its position manually!" << std::endl;
202  }
203  prefix_ = r.is_prefix;
204  decoy_string_ = r.name;
205  // decoy string and position was extracted successfully
206  OPENMS_LOG_INFO << "Using " << (prefix_ ? "prefix" : "suffix") << " decoy string '" << decoy_string_ << "'" << std::endl;
207  }
208 
209  //---------------------------------------------------------------
210  // parsing parameters, correcting xtandem and MSGFPlus parameters
211  //---------------------------------------------------------------
212  ProteaseDigestion enzyme;
213  enzyme.setEnzyme(enzyme_name_);
214  enzyme.setSpecificity(enzyme.getSpecificityByName(enzyme_specificity_));
215 
216  bool xtandem_fix_parameters = true;
217  bool msgfplus_fix_parameters = true;
218 
219  // determine if search engine is solely xtandem or MSGFPlus
220  for (const auto& prot_id : prot_ids)
221  {
222  String search_engine = prot_id.getSearchEngine();
223  StringUtils::toUpper(search_engine);
224  if (search_engine != "XTANDEM") { xtandem_fix_parameters = false; }
225  if (!(search_engine == "MSGFPLUS" || search_engine == "MS-GF+")) { msgfplus_fix_parameters = false; }
226  }
227 
228  // solely MSGFPlus -> Trypsin/P as enzyme
229  if (msgfplus_fix_parameters && enzyme.getEnzymeName() == "Trypsin")
230  {
231  OPENMS_LOG_WARN << "MSGFPlus detected but enzyme cutting rules were set to Trypsin. Correcting to Trypsin/P to copy with special cutting rule in MSGFPlus." << std::endl;
232  enzyme.setEnzyme("Trypsin/P");
233  }
234 
235  //-------------------------------------------------------------
236  // calculations
237  //-------------------------------------------------------------
238  // cache the first proteins
239  const size_t PROTEIN_CACHE_SIZE = 4e5; // 400k should be enough for most DB's and is not too hard on memory either (~200 MB FASTA)
240 
241  this->startProgress(0, 1, "Load first DB chunk");
242  proteins.cacheChunk(PROTEIN_CACHE_SIZE);
243  this->endProgress();
244 
245  if (proteins.empty()) // we do not allow an empty database
246  {
247  OPENMS_LOG_ERROR << "Error: An empty database was provided. Mapping makes no sense. Aborting..." << std::endl;
248  return DATABASE_EMPTY;
249  }
250 
251  if (pep_ids.empty()) // Aho-Corasick requires non-empty input; but we allow this case, since the TOPP tool should not crash when encountering a bad raw file (with no PSMs)
252  {
253  OPENMS_LOG_WARN << "Warning: An empty set of peptide identifications was provided. Output will be empty as well." << std::endl;
254  if (!keep_unreferenced_proteins_)
255  {
256  // delete only protein hits, not whole ID runs incl. meta data:
257  for (std::vector<ProteinIdentification>::iterator it = prot_ids.begin();
258  it != prot_ids.end(); ++it)
259  {
260  it->getHits().clear();
261  }
262  }
263  return PEPTIDE_IDS_EMPTY;
264  }
265 
266  FoundProteinFunctor func(enzyme, xtandem_fix_parameters); // store the matches
267  Map<String, Size> acc_to_prot; // map: accessions --> FASTA protein index
268  std::vector<bool> protein_is_decoy; // protein index -> is decoy?
269  std::vector<std::string> protein_accessions; // protein index -> accession
270 
271  bool invalid_protein_sequence = false; // check for proteins with modifications, i.e. '[' or '(', and throw an exception
272 
273  { // new scope - forget data after search
274 
275  /*
276  BUILD Peptide DB
277  */
278  bool has_illegal_AAs(false);
280  for (std::vector<PeptideIdentification>::const_iterator it1 = pep_ids.begin(); it1 != pep_ids.end(); ++it1)
281  {
282  //String run_id = it1->getIdentifier();
283  const std::vector<PeptideHit>& hits = it1->getHits();
284  for (std::vector<PeptideHit>::const_iterator it2 = hits.begin(); it2 != hits.end(); ++it2)
285  {
286  //
287  // Warning:
288  // do not skip over peptides here, since the results are iterated in the same way
289  //
290  String seq = it2->getSequence().toUnmodifiedString().remove('*'); // make a copy, i.e. do NOT change the peptide sequence!
291  if (seqan::isAmbiguous(seqan::AAString(seq.c_str())))
292  { // do not quit here, to show the user all sequences .. only quit after loop
293  OPENMS_LOG_ERROR << "Peptide sequence '" << it2->getSequence() << "' contains one or more ambiguous amino acids (B|J|Z|X).\n";
294  has_illegal_AAs = true;
295  }
296  if (IL_equivalent_) // convert L to I;
297  {
298  seq.substitute('L', 'I');
299  }
300  appendValue(pep_DB, seq.c_str());
301  }
302  }
303  if (has_illegal_AAs)
304  {
305  OPENMS_LOG_ERROR << "One or more peptides contained illegal amino acids. This is not allowed!"
306  << "\nPlease either remove the peptide or replace it with one of the unambiguous ones (while allowing for ambiguous AA's to match the protein)." << std::endl;;
307  }
308 
309  OPENMS_LOG_INFO << "Mapping " << length(pep_DB) << " peptides to " << (proteins.size() == PROTEIN_CACHE_SIZE ? "? (unknown number of)" : String(proteins.size())) << " proteins." << std::endl;
310 
311  if (length(pep_DB) == 0)
312  { // Aho-Corasick will crash if given empty needles as input
313  OPENMS_LOG_WARN << "Warning: Peptide identifications have no hits inside! Output will be empty as well." << std::endl;
314  return PEPTIDE_IDS_EMPTY;
315  }
316 
317  /*
318  Aho Corasick (fast)
319  */
320  OPENMS_LOG_INFO << "Searching with up to " << aaa_max_ << " ambiguous amino acid(s) and " << mm_max_ << " mismatch(es)!" << std::endl;
322  OPENMS_LOG_INFO << "Building trie ...";
323  StopWatch s;
324  s.start();
326  AhoCorasickAmbiguous::initPattern(pep_DB, aaa_max_, mm_max_, pattern);
327  s.stop();
328  OPENMS_LOG_INFO << " done (" << int(s.getClockTime()) << "s)" << std::endl;
329  s.reset();
330 
331  uint16_t count_j_proteins(0);
332  bool has_active_data = true; // becomes false if end of FASTA file is reached
333  const std::string jumpX(aaa_max_ + mm_max_ + 1, 'X'); // jump over stretches of 'X' which cost a lot of time; +1 because AXXA is a valid hit for aaa_max == 2 (cannot split it)
334  // use very large target value for progress if DB size is unknown (did not fit into first chunk)
335  this->startProgress(0, proteins.size() == PROTEIN_CACHE_SIZE ? std::numeric_limits<SignedSize>::max() : proteins.size(), "Aho-Corasick");
336  std::atomic<int> progress_prots(0);
337 #ifdef _OPENMP
338 #pragma omp parallel
339 #endif
340  {
341  FoundProteinFunctor func_threads(enzyme, xtandem_fix_parameters);
342  Map<String, Size> acc_to_prot_thread; // map: accessions --> FASTA protein index
343  AhoCorasickAmbiguous fuzzyAC;
344  String prot;
345 
346  while (true)
347  {
348  #pragma omp barrier // all threads need to be here, since we are about to swap protein data
349  #pragma omp single
350  {
351  DEBUG_ONLY std::cerr << " activating cache ...\n";
352  has_active_data = proteins.activateCache(); // swap in last cache
353  protein_accessions.resize(proteins.getChunkOffset() + proteins.chunkSize());
354  } // implicit barrier here
355 
356  if (!has_active_data) break; // leave while-loop
357  SignedSize prot_count = (SignedSize)proteins.chunkSize();
358 
359  #pragma omp master
360  {
361  DEBUG_ONLY std::cerr << "Filling Protein Cache ...";
362  proteins.cacheChunk(PROTEIN_CACHE_SIZE);
363  protein_is_decoy.resize(proteins.getChunkOffset() + prot_count);
364  for (SignedSize i = 0; i < prot_count; ++i)
365  { // do this in master only, to avoid false sharing
366  const String& seq = proteins.chunkAt(i).identifier;
367  protein_is_decoy[i + proteins.getChunkOffset()] = (prefix_ ? seq.hasPrefix(decoy_string_) : seq.hasSuffix(decoy_string_));
368  }
369  DEBUG_ONLY std::cerr << " done" << std::endl;
370  }
371  DEBUG_ONLY std::cerr << " starting for loop \n";
372  // search all peptides in each protein
373  #pragma omp for schedule(dynamic, 100) nowait
374  for (SignedSize i = 0; i < prot_count; ++i)
375  {
376  ++progress_prots; // atomic
377  if (omp_get_thread_num() == 0)
378  {
379  this->setProgress(progress_prots);
380  }
381 
382  prot = proteins.chunkAt(i).sequence;
383  prot.remove('*');
384 
385  // check for invalid sequences with modifications
386  if (prot.has('[') || prot.has('('))
387  {
388  invalid_protein_sequence = true; // not omp-critical because its write-only
389  // we cannot throw an exception here, since we'd need to catch it within the parallel region
390  }
391 
392  // convert L/J to I; also replace 'J' in proteins
393  if (IL_equivalent_)
394  {
395  prot.substitute('L', 'I');
396  prot.substitute('J', 'I');
397  }
398  else
399  { // warn if 'J' is found (it eats into aaa_max)
400  if (prot.has('J'))
401  {
402  #pragma omp atomic
403  ++count_j_proteins;
404  }
405  }
406 
407  Size prot_idx = i + proteins.getChunkOffset();
408 
409  // test if protein was a hit
410  Size hits_total = func_threads.filter_passed + func_threads.filter_rejected;
411 
412  // check if there are stretches of 'X'
413  if (prot.has('X'))
414  {
415  // create chunks of the protein (splitting it at stretches of 'X..X') and feed them to AC one by one
416  size_t offset = -1, start = 0;
417  while ((offset = prot.find(jumpX, offset + 1)) != std::string::npos)
418  {
419  //std::cout << "found X..X at " << offset << " in protein " << proteins[i].identifier << "\n";
420  addHits_(fuzzyAC, pattern, pep_DB, prot.substr(start, offset + jumpX.size() - start), prot, prot_idx, (int)start, func_threads);
421  // skip ahead while we encounter more X...
422  while (offset + jumpX.size() < prot.size() && prot[offset + jumpX.size()] == 'X') ++offset;
423  start = offset;
424  //std::cout << " new start: " << start << "\n";
425  }
426  // last chunk
427  if (start < prot.size())
428  {
429  addHits_(fuzzyAC, pattern, pep_DB, prot.substr(start), prot, prot_idx, (int)start, func_threads);
430  }
431  }
432  else
433  {
434  addHits_(fuzzyAC, pattern, pep_DB, prot, prot, prot_idx, 0, func_threads);
435  }
436  // was protein found?
437  if (hits_total < func_threads.filter_passed + func_threads.filter_rejected)
438  {
439  protein_accessions[prot_idx] = proteins.chunkAt(i).identifier;
440  acc_to_prot_thread[protein_accessions[prot_idx]] = prot_idx;
441  }
442  } // end parallel FOR
443 
444  // join results again
445  DEBUG_ONLY std::cerr << " critical now \n";
446  #ifdef _OPENMP
447  #pragma omp critical(PeptideIndexer_joinAC)
448  #endif
449  {
450  s.start();
451  // hits
452  func.merge(func_threads);
453  // accession -> index
454  acc_to_prot.insert(acc_to_prot_thread.begin(), acc_to_prot_thread.end());
455  acc_to_prot_thread.clear();
456  s.stop();
457  } // OMP end critical
458  } // end readChunk
459  } // OMP end parallel
460  this->endProgress();
461  std::cout << "Merge took: " << s.toString() << "\n";
462  mu.after();
463  std::cout << mu.delta("Aho-Corasick") << "\n\n";
464 
465  OPENMS_LOG_INFO << "\nAho-Corasick done:\n found " << func.filter_passed << " hits for " << func.pep_to_prot.size() << " of " << length(pep_DB) << " peptides.\n";
466 
467  // write some stats
468  OPENMS_LOG_INFO << "Peptide hits passing enzyme filter: " << func.filter_passed << "\n"
469  << " ... rejected by enzyme filter: " << func.filter_rejected << std::endl;
470 
471  if (count_j_proteins)
472  {
473  OPENMS_LOG_WARN << "PeptideIndexer found " << count_j_proteins << " protein sequences in your database containing the amino acid 'J'."
474  << "To match 'J' in a protein, an ambiguous amino acid placeholder for I/L will be used.\n"
475  << "This costs runtime and eats into the 'aaa_max' limit, leaving less opportunity for B/Z/X matches.\n"
476  << "If you want 'J' to be treated as unambiguous, enable '-IL_equivalent'!" << std::endl;
477  }
478 
479  } // end local scope
480 
481  //
482  // do mapping
483  //
484  // index existing proteins
485  Map<String, Size> runid_to_runidx; // identifier to index
486  for (Size run_idx = 0; run_idx < prot_ids.size(); ++run_idx)
487  {
488  runid_to_runidx[prot_ids[run_idx].getIdentifier()] = run_idx;
489  }
490 
491  // for peptides --> proteins
492  Size stats_matched_unique(0);
493  Size stats_matched_multi(0);
494  Size stats_unmatched(0); // no match to DB
495  Size stats_count_m_t(0); // match to Target DB
496  Size stats_count_m_d(0); // match to Decoy DB
497  Size stats_count_m_td(0); // match to T+D DB
498 
499  Map<Size, std::set<Size> > runidx_to_protidx; // in which protID do appear which proteins (according to mapped peptides)
500 
501  Size pep_idx(0);
502  for (std::vector<PeptideIdentification>::iterator it1 = pep_ids.begin(); it1 != pep_ids.end(); ++it1)
503  {
504  // which ProteinIdentification does the peptide belong to?
505  Size run_idx = runid_to_runidx[it1->getIdentifier()];
506 
507  std::vector<PeptideHit>& hits = it1->getHits();
508 
509  for (std::vector<PeptideHit>::iterator it2 = hits.begin(); it2 != hits.end(); ++it2)
510  {
511  // clear protein accessions
512  it2->setPeptideEvidences(std::vector<PeptideEvidence>());
513 
514  //
515  // is this a decoy hit?
516  //
517  bool matches_target(false);
518  bool matches_decoy(false);
519 
520  std::set<Size> prot_indices;
521  // add new protein references
522  for (std::set<PeptideProteinMatchInformation>::const_iterator it_i = func.pep_to_prot[pep_idx].begin();
523  it_i != func.pep_to_prot[pep_idx].end(); ++it_i)
524  {
525  prot_indices.insert(it_i->protein_index);
526  const String& accession = protein_accessions[it_i->protein_index];
527  PeptideEvidence pe(accession, it_i->position, it_i->position + (int)it2->getSequence().size() - 1, it_i->AABefore, it_i->AAAfter);
528  it2->addPeptideEvidence(pe);
529 
530  runidx_to_protidx[run_idx].insert(it_i->protein_index); // fill protein hits
531 
532  if (protein_is_decoy[it_i->protein_index])
533  {
534  matches_decoy = true;
535  }
536  else
537  {
538  matches_target = true;
539  }
540  }
541 
542  if (matches_decoy && matches_target)
543  {
544  it2->setMetaValue("target_decoy", "target+decoy");
545  ++stats_count_m_td;
546  }
547  else if (matches_target)
548  {
549  it2->setMetaValue("target_decoy", "target");
550  ++stats_count_m_t;
551  }
552  else if (matches_decoy)
553  {
554  it2->setMetaValue("target_decoy", "decoy");
555  ++stats_count_m_d;
556  } // else: could match to no protein (i.e. both are false)
557  //else ... // not required (handled below; see stats_unmatched);
558 
559  if (prot_indices.size() == 1)
560  {
561  it2->setMetaValue("protein_references", "unique");
562  ++stats_matched_unique;
563  }
564  else if (prot_indices.size() > 1)
565  {
566  it2->setMetaValue("protein_references", "non-unique");
567  ++stats_matched_multi;
568  }
569  else
570  {
571  it2->setMetaValue("protein_references", "unmatched");
572  ++stats_unmatched;
573  if (stats_unmatched < 15) OPENMS_LOG_INFO << "Unmatched peptide: " << it2->getSequence() << "\n";
574  else if (stats_unmatched == 15) OPENMS_LOG_INFO << "Unmatched peptide: ...\n";
575  }
576 
577  ++pep_idx; // next hit
578  }
579 
580  }
581 
582  Size total_peptides = stats_count_m_t + stats_count_m_d + stats_count_m_td + stats_unmatched;
583  OPENMS_LOG_INFO << "-----------------------------------\n";
584  OPENMS_LOG_INFO << "Peptide statistics\n";
585  OPENMS_LOG_INFO << "\n";
586  OPENMS_LOG_INFO << " unmatched : " << stats_unmatched << " (" << stats_unmatched * 100 / total_peptides << " %)\n";
587  OPENMS_LOG_INFO << " target/decoy:\n";
588  OPENMS_LOG_INFO << " match to target DB only: " << stats_count_m_t << " (" << stats_count_m_t * 100 / total_peptides << " %)\n";
589  OPENMS_LOG_INFO << " match to decoy DB only : " << stats_count_m_d << " (" << stats_count_m_d * 100 / total_peptides << " %)\n";
590  OPENMS_LOG_INFO << " match to both : " << stats_count_m_td << " (" << stats_count_m_td * 100 / total_peptides << " %)\n";
591  OPENMS_LOG_INFO << "\n";
592  OPENMS_LOG_INFO << " mapping to proteins:\n";
593  OPENMS_LOG_INFO << " no match (to 0 protein) : " << stats_unmatched << "\n";
594  OPENMS_LOG_INFO << " unique match (to 1 protein) : " << stats_matched_unique << "\n";
595  OPENMS_LOG_INFO << " non-unique match (to >1 protein): " << stats_matched_multi << std::endl;
596 
598  Size stats_matched_proteins(0), stats_matched_new_proteins(0), stats_orphaned_proteins(0), stats_proteins_target(0), stats_proteins_decoy(0);
599 
600  // all peptides contain the correct protein hit references, now update the protein hits
601  for (Size run_idx = 0; run_idx < prot_ids.size(); ++run_idx)
602  {
603  std::set<Size> masterset = runidx_to_protidx[run_idx]; // all protein matches from above
604 
605  std::vector<ProteinHit>& phits = prot_ids[run_idx].getHits();
606  {
607  // go through existing protein hits and count orphaned proteins (with no peptide hits)
608  std::vector<ProteinHit> orphaned_hits;
609  for (std::vector<ProteinHit>::iterator p_hit = phits.begin(); p_hit != phits.end(); ++p_hit)
610  {
611  const String& acc = p_hit->getAccession();
612  if (!acc_to_prot.has(acc)) // acc_to_prot only contains found proteins from current run
613  { // old hit is orphaned
614  ++stats_orphaned_proteins;
615  if (keep_unreferenced_proteins_)
616  {
617  p_hit->setMetaValue("target_decoy", "");
618  orphaned_hits.push_back(*p_hit);
619  }
620  }
621  }
622  // only keep orphaned hits (if any)
623  phits = orphaned_hits;
624  }
625 
626  // add new protein hits
628  phits.reserve(phits.size() + masterset.size());
629  for (std::set<Size>::const_iterator it = masterset.begin(); it != masterset.end(); ++it)
630  {
631  ProteinHit hit;
632  hit.setAccession(protein_accessions[*it]);
633 
634  if (write_protein_sequence_ || write_protein_description_)
635  {
636  proteins.readAt(fe, *it);
637  if (write_protein_sequence_)
638  {
639  hit.setSequence(fe.sequence);
640  } // no else, since sequence is empty by default
641  if (write_protein_description_)
642  {
643  hit.setDescription(fe.description);
644  } // no else, since description is empty by default
645  }
646  if (protein_is_decoy[*it])
647  {
648  hit.setMetaValue("target_decoy", "decoy");
649  ++stats_proteins_decoy;
650  }
651  else
652  {
653  hit.setMetaValue("target_decoy", "target");
654  ++stats_proteins_target;
655  }
656  phits.push_back(hit);
657  ++stats_matched_new_proteins;
658  }
659  stats_matched_proteins += phits.size();
660  }
661 
662 
663  OPENMS_LOG_INFO << "-----------------------------------\n";
664  OPENMS_LOG_INFO << "Protein statistics\n";
665  OPENMS_LOG_INFO << "\n";
666  OPENMS_LOG_INFO << " total proteins searched: " << proteins.size() << "\n";
667  OPENMS_LOG_INFO << " matched proteins : " << stats_matched_proteins << " (" << stats_matched_new_proteins << " new)\n";
668  if (stats_matched_proteins)
669  { // prevent Division-by-0 Exception
670  OPENMS_LOG_INFO << " matched target proteins: " << stats_proteins_target << " (" << stats_proteins_target * 100 / stats_matched_proteins << " %)\n";
671  OPENMS_LOG_INFO << " matched decoy proteins : " << stats_proteins_decoy << " (" << stats_proteins_decoy * 100 / stats_matched_proteins << " %)\n";
672  }
673  OPENMS_LOG_INFO << " orphaned proteins : " << stats_orphaned_proteins << (keep_unreferenced_proteins_ ? " (all kept)" : " (all removed)\n");
674  OPENMS_LOG_INFO << "-----------------------------------" << std::endl;
675 
676 
678  bool has_error = false;
679 
680  if (invalid_protein_sequence)
681  {
682  OPENMS_LOG_ERROR << "Error: One or more protein sequences contained the characters '[' or '(', which are illegal in protein sequences."
683  << "\nPeptide hits might be masked by these characters (which usually indicate presence of modifications).\n";
684  has_error = true;
685  }
686 
687  if ((stats_count_m_d + stats_count_m_td) == 0)
688  {
689  String msg("No peptides were matched to the decoy portion of the database! Did you provide the correct concatenated database? Are your 'decoy_string' (=" + String(decoy_string_) + ") and 'decoy_string_position' (=" + String(param_.getValue("decoy_string_position")) + ") settings correct?");
690  if (missing_decoy_action_ == "error")
691  {
692  OPENMS_LOG_ERROR << "Error: " << msg << "\nSet 'missing_decoy_action' to 'warn' if you are sure this is ok!\nAborting ..." << std::endl;
693  has_error = true;
694  }
695  else if (missing_decoy_action_ == "warn")
696  {
697  OPENMS_LOG_WARN << "Warn: " << msg << "\nSet 'missing_decoy_action' to 'error' if you want to elevate this to an error!" << std::endl;
698  }
699  else // silent
700  {
701  }
702  }
703 
704  if ((!allow_unmatched_) && (stats_unmatched > 0))
705  {
706  OPENMS_LOG_ERROR << "PeptideIndexer found unmatched peptides, which could not be associated to a protein.\n"
707  << "Potential solutions:\n"
708  << " - check your FASTA database for completeness\n"
709  << " - set 'enzyme:specificity' to match the identification parameters of the search engine\n"
710  << " - some engines (e.g. X! Tandem) employ loose cutting rules generating non-tryptic peptides;\n"
711  << " if you trust them, disable enzyme specificity\n"
712  << " - increase 'aaa_max' to allow more ambiguous amino acids\n"
713  << " - as a last resort: use the 'allow_unmatched' option to accept unmatched peptides\n"
714  << " (note that unmatched peptides cannot be used for FDR calculation or quantification)\n";
715  has_error = true;
716  }
717 
718  if (has_error)
719  {
720  OPENMS_LOG_ERROR << "Result files will be written, but PeptideIndexer will exit with an error code." << std::endl;
721  return UNEXPECTED_RESULT;
722  }
723  return EXECUTION_OK;
724  }
725 
726  const String& getDecoyString() const;
727 
728  bool isPrefix() const;
729 
730  protected:
731 
733  {
734  OpenMS::Size protein_index; //< index of the protein the peptide is contained in
735  OpenMS::Int position; //< the position of the peptide in the protein
736  char AABefore; //< the amino acid after the peptide in the protein
737  char AAAfter; //< the amino acid before the peptide in the protein
738 
739  const std::tuple<const Size&, const Int&, const char&, const char&> tie() const
740  {
741  return std::tie(protein_index, position, AABefore, AAAfter);
742  }
743  bool operator<(const PeptideProteinMatchInformation& other) const
744  {
745  return tie() < other.tie();
746  }
748  {
749  return tie() == other.tie();
750  }
751  };
752 
754  {
755  public:
756  typedef std::map<OpenMS::Size, std::set<PeptideProteinMatchInformation> > MapType;
757  MapType pep_to_prot; //< peptide index --> protein indices
758  OpenMS::Size filter_passed; //< number of accepted hits (passing addHit() constraints)
759  OpenMS::Size filter_rejected; //< number of rejected hits (not passing addHit())
760 
761  private:
763  bool xtandem_; //< are we checking xtandem cleavage rules?
764 
765  public:
766  explicit FoundProteinFunctor(const ProteaseDigestion& enzyme, bool xtandem) :
767  pep_to_prot(), filter_passed(0), filter_rejected(0), enzyme_(enzyme), xtandem_(xtandem)
768  {
769  }
770 
772  {
773  if (pep_to_prot.empty())
774  { // first merge is easy
775  pep_to_prot.swap(other.pep_to_prot);
776  }
777  else
778  {
779  for (FoundProteinFunctor::MapType::const_iterator it = other.pep_to_prot.begin(); it != other.pep_to_prot.end(); ++it)
780  { // augment set
781  this->pep_to_prot[it->first].insert(other.pep_to_prot[it->first].begin(), other.pep_to_prot[it->first].end());
782  }
783  other.pep_to_prot.clear();
784  }
785  // cheap members
786  this->filter_passed += other.filter_passed;
787  other.filter_passed = 0;
788  this->filter_rejected += other.filter_rejected;
789  other.filter_rejected = 0;
790  }
791 
792  void addHit(const OpenMS::Size idx_pep,
793  const OpenMS::Size idx_prot,
794  const OpenMS::Size len_pep,
795  const OpenMS::String& seq_prot,
797  {
798  if (enzyme_.isValidProduct(seq_prot, position, len_pep, true, true, xtandem_))
799  {
801  match.protein_index = idx_prot;
802  match.position = position;
803  match.AABefore = (position == 0) ? PeptideEvidence::N_TERMINAL_AA : seq_prot[position - 1];
804  match.AAAfter = (position + len_pep >= seq_prot.size()) ? PeptideEvidence::C_TERMINAL_AA : seq_prot[position + len_pep];
805  pep_to_prot[idx_pep].insert(match);
806  ++filter_passed;
807  }
808  else
809  {
810  //std::cerr << "REJECTED Peptide " << seq_pep << " with hit to protein "
811  // << seq_prot << " at position " << position << std::endl;
812  ++filter_rejected;
813  }
814  }
815 
816  };
817 
818  inline void addHits_(AhoCorasickAmbiguous& fuzzyAC, const AhoCorasickAmbiguous::FuzzyACPattern& pattern, const AhoCorasickAmbiguous::PeptideDB& pep_DB, const String& prot, const String& full_prot, SignedSize idx_prot, Int offset, FoundProteinFunctor& func_threads) const
819  {
820  fuzzyAC.setProtein(prot);
821  while (fuzzyAC.findNext(pattern))
822  {
823  const seqan::Peptide& tmp_pep = pep_DB[fuzzyAC.getHitDBIndex()];
824  func_threads.addHit(fuzzyAC.getHitDBIndex(), idx_prot, length(tmp_pep), full_prot, fuzzyAC.getHitProteinPosition() + offset);
825  }
826  }
827 
828  void updateMembers_() override;
829 
831  bool prefix_;
835 
841 
844  };
845 }
846 
LogStream.h
DefaultParamHandler.h
OpenMS::TOPPBase
Base class for TOPP applications.
Definition: TOPPBase.h:144
OpenMS::Param::copy
Param copy(const String &prefix, bool remove_prefix=false) const
Returns a new Param object containing all entries that start with prefix.
FileHandler.h
OpenMS::PeptideIndexing::addHits_
void addHits_(AhoCorasickAmbiguous &fuzzyAC, const AhoCorasickAmbiguous::FuzzyACPattern &pattern, const AhoCorasickAmbiguous::PeptideDB &pep_DB, const String &prot, const String &full_prot, SignedSize idx_prot, Int offset, FoundProteinFunctor &func_threads) const
Definition: PeptideIndexing.h:818
OpenMS::PeptideIndexing::PEPTIDE_IDS_EMPTY
Definition: PeptideIndexing.h:134
OpenMS::PeptideIndexing::enzyme_name_
String enzyme_name_
Definition: PeptideIndexing.h:833
OpenMS::SysInfo::MemUsage
A convenience class to report either absolute or delta (between two timepoints) RAM usage.
Definition: SysInfo.h:83
OpenMS::ProteinHit::setDescription
void setDescription(const String &description)
sets the description of the protein
OpenMS::StopWatch::stop
bool stop()
OpenMS::IdXMLFile::store
void store(const String &filename, const std::vector< ProteinIdentification > &protein_ids, const std::vector< PeptideIdentification > &peptide_ids, const String &document_id="")
Stores the data in an idXML file.
OpenMS::ProteinHit::setSequence
void setSequence(const String &sequence)
sets the protein sequence
OpenMS::ProteaseDigestion::setEnzyme
void setEnzyme(const String &name)
Sets the enzyme for the digestion (by name)
OpenMS::PeptideIndexing::FoundProteinFunctor::filter_passed
OpenMS::Size filter_passed
Definition: PeptideIndexing.h:758
OpenMS::AhoCorasickAmbiguous::findNext
bool findNext(const FuzzyACPattern &pattern)
Enumerate hits.
Definition: AhoCorasickAmbiguous.h:1037
StopWatch.h
OpenMS::AhoCorasickAmbiguous::setProtein
void setProtein(const String &protein_sequence)
Reset to new protein sequence. All previous data is forgotten.
Definition: AhoCorasickAmbiguous.h:1024
OpenMS::PeptideIndexing::FoundProteinFunctor::enzyme_
ProteaseDigestion enzyme_
Definition: PeptideIndexing.h:762
OpenMS::String::substitute
String & substitute(char from, char to)
Replaces all occurrences of the character from by the character to.
OpenMS::Map::has
bool has(const Key &key) const
Test whether the map contains the given key.
Definition: Map.h:108
OpenMS::FASTAFile::FASTAEntry::sequence
String sequence
Definition: FASTAFile.h:80
OpenMS::PeptideIndexing::decoy_string_
String decoy_string_
Definition: PeptideIndexing.h:830
OpenMS::String
A more convenient string class.
Definition: String.h:58
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::operator==
bool operator==(const PeptideProteinMatchInformation &other) const
Definition: PeptideIndexing.h:747
PeptideEvidence.h
FASTAContainer.h
OpenMS::PeptideIndexing::FoundProteinFunctor::addHit
void addHit(const OpenMS::Size idx_pep, const OpenMS::Size idx_prot, const OpenMS::Size len_pep, const OpenMS::String &seq_prot, OpenMS::Int position)
Definition: PeptideIndexing.h:792
OpenMS::FASTAFile::FASTAEntry::description
String description
Definition: FASTAFile.h:79
OpenMS::Size
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
IdXMLFile.h
OpenMS::ProteinHit
Representation of a protein hit.
Definition: ProteinHit.h:57
OpenMS::AhoCorasickAmbiguous::getHitProteinPosition
Int getHitProteinPosition()
Offset into protein sequence where hit was found.
Definition: AhoCorasickAmbiguous.h:1057
seqan::AAString
String< AAcid, Alloc< void > > AAString
Definition: AhoCorasickAmbiguous.h:206
OpenMS::AhoCorasickAmbiguous
Extended Aho-Corasick algorithm capable of matching ambiguous amino acids in the pattern (i....
Definition: AhoCorasickAmbiguous.h:970
OpenMS::PeptideIndexing::PeptideProteinMatchInformation
Definition: PeptideIndexing.h:732
OpenMS::Param::getValue
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
OpenMS::PeptideIndexing::run
ExitCodes run(std::vector< FASTAFile::FASTAEntry > &proteins, std::vector< ProteinIdentification > &prot_ids, std::vector< PeptideIdentification > &pep_ids)
forward for old interface and pyOpenMS; use run<T>() for more control
Definition: PeptideIndexing.h:147
OpenMS::SysInfo::MemUsage::delta
String delta(const String &event="delta")
DEBUG_ONLY
#define DEBUG_ONLY
Definition: AhoCorasickAmbiguous.h:46
OpenMS::AhoCorasickAmbiguous::initPattern
static void initPattern(const PeptideDB &pep_db, const int aaa_max, const int mm_max, FuzzyACPattern &pattern)
Construct a trie from a set of peptide sequences (which are to be found in a protein).
Definition: AhoCorasickAmbiguous.h:991
OPENMS_LOG_WARN
#define OPENMS_LOG_WARN
Macro if a warning, a piece of information which should be read by the user, should be logged.
Definition: LogStream.h:460
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::protein_index
OpenMS::Size protein_index
Definition: PeptideIndexing.h:734
OpenMS::EnzymaticDigestion::getSpecificityByName
static Specificity getSpecificityByName(const String &name)
OpenMS::PeptideIndexing::FoundProteinFunctor::filter_rejected
OpenMS::Size filter_rejected
Definition: PeptideIndexing.h:759
OpenMS::FASTAFile::FASTAEntry
FASTA entry type (identifier, description and sequence)
Definition: FASTAFile.h:76
OpenMS::PeptideIndexing::UNEXPECTED_RESULT
Definition: PeptideIndexing.h:136
OpenMS::String::hasPrefix
bool hasPrefix(const String &string) const
true if String begins with string, false otherwise
OpenMS::PeptideIndexing::FoundProteinFunctor::pep_to_prot
MapType pep_to_prot
Definition: PeptideIndexing.h:757
ListUtils.h
OpenMS::DefaultParamHandler
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:91
OpenMS::IdXMLFile::load
void load(const String &filename, std::vector< ProteinIdentification > &protein_ids, std::vector< PeptideIdentification > &peptide_ids)
Loads the identifications of an idXML file without identifier.
OpenMS::ProteinHit::setAccession
void setAccession(const String &accession)
sets the accession of the protein
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::operator<
bool operator<(const PeptideProteinMatchInformation &other) const
Definition: PeptideIndexing.h:743
OpenMS::PeptideEvidence::N_TERMINAL_AA
static const char N_TERMINAL_AA
Definition: PeptideEvidence.h:60
OpenMS::PeptideIndexing::run
ExitCodes run(FASTAContainer< T > &proteins, std::vector< ProteinIdentification > &prot_ids, std::vector< PeptideIdentification > &pep_ids)
Re-index peptide identifications honoring enzyme cutting rules, ambiguous amino acids and target/deco...
Definition: PeptideIndexing.h:189
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
OpenMS::DataValue::toBool
bool toBool() const
Conversion to bool.
OpenMS::AhoCorasickAmbiguous::getHitDBIndex
Size getHitDBIndex()
Get index of hit into peptide database of the pattern.
Definition: AhoCorasickAmbiguous.h:1047
ProteaseDigestion.h
OpenMS::PeptideIndexing::FoundProteinFunctor::merge
void merge(FoundProteinFunctor &other)
Definition: PeptideIndexing.h:771
OpenMS::PeptideIndexing::FoundProteinFunctor::FoundProteinFunctor
FoundProteinFunctor(const ProteaseDigestion &enzyme, bool xtandem)
Definition: PeptideIndexing.h:766
OpenMS::ProgressLogger
Base class for all classes that want to report their progress.
Definition: ProgressLogger.h:54
ProteaseDB.h
OpenMS::MetaInfoInterface::setMetaValue
void setMetaValue(const String &name, const DataValue &value)
Sets the DataValue corresponding to a name.
ProgressLogger.h
FASTAFile.h
OpenMS::FASTAContainer< TFI_Vector >
FASTAContainer<TFI_Vector> simply takes an existing vector of FASTAEntries and provides the same inte...
Definition: FASTAContainer.h:243
OpenMS::StopWatch::getClockTime
double getClockTime() const
OpenMS::StopWatch::reset
void reset()
int
OpenMS::PeptideIndexing::aaa_max_
Int aaa_max_
Definition: PeptideIndexing.h:842
OpenMS::String::has
bool has(Byte byte) const
true if String contains the byte, false otherwise
ProteinIdentification.h
OpenMS::PeptideEvidence
Representation of a peptide evidence.
Definition: PeptideEvidence.h:50
OpenMS::PeptideIndexing::write_protein_sequence_
bool write_protein_sequence_
Definition: PeptideIndexing.h:836
OpenMS::AhoCorasickAmbiguous::FuzzyACPattern
::seqan::Pattern< PeptideDB, ::seqan::FuzzyAC > FuzzyACPattern
Definition: AhoCorasickAmbiguous.h:974
OpenMS::DefaultParamHandler::setParameters
void setParameters(const Param &param)
Sets the parameters.
OpenMS::String::hasSuffix
bool hasSuffix(const String &string) const
true if String ends with string, false otherwise
OpenMS::PeptideIndexing::missing_decoy_action_
String missing_decoy_action_
Definition: PeptideIndexing.h:832
AhoCorasickAmbiguous.h
OpenMS::DefaultParamHandler::getParameters
const Param & getParameters() const
Non-mutable access to the parameters.
OpenMS::AhoCorasickAmbiguous::PeptideDB
::seqan::StringSet<::seqan::AAString > PeptideDB
Definition: AhoCorasickAmbiguous.h:973
OpenMS::StopWatch::toString
String toString() const
get a compact representation of the current time status.
OpenMS::PeptideIndexing::ExitCodes
ExitCodes
Exit codes.
Definition: PeptideIndexing.h:130
SysInfo.h
OpenMS::SysInfo::MemUsage::after
void after()
record data for the second timepoint
OpenMS::StopWatch
StopWatch Class.
Definition: StopWatch.h:59
OpenMS::PeptideIndexing::prefix_
bool prefix_
Definition: PeptideIndexing.h:831
OpenMS::ProteaseDigestion::isValidProduct
bool isValidProduct(const String &protein, int pep_pos, int pep_length, bool ignore_missed_cleavages=true, bool allow_nterm_protein_cleavage=false, bool allow_random_asp_pro_cleavage=false) const
Variant of EnzymaticDigestion::isValidProduct() with support for n-term protein cleavage and random D...
StringUtils.h
main
int main(int argc, const char **argv)
Definition: INIFileEditor.cpp:73
OpenMS::PeptideIndexing::FoundProteinFunctor::MapType
std::map< OpenMS::Size, std::set< PeptideProteinMatchInformation > > MapType
Definition: PeptideIndexing.h:756
OpenMS::StopWatch::start
bool start()
OpenMS::PeptideIndexing::write_protein_description_
bool write_protein_description_
Definition: PeptideIndexing.h:837
OpenMS::Param::update
bool update(const Param &p_outdated, const bool add_unknown=false)
Rescue parameter values from p_outdated to current param.
OpenMS::SignedSize
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
OpenMS::EnzymaticDigestion::getEnzymeName
String getEnzymeName() const
Returns the enzyme for the digestion.
OPENMS_LOG_ERROR
#define OPENMS_LOG_ERROR
Macro to be used if non-fatal error are reported (processing continues)
Definition: LogStream.h:455
OpenMS::OpenMS_Log_debug
Logger::LogStream OpenMS_Log_debug
Global static instance of a LogStream to capture messages classified as debug output....
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::AAAfter
char AAAfter
Definition: PeptideIndexing.h:737
OpenMS::String::substr
String substr(size_t pos=0, size_t n=npos) const
Wrapper for the STL substr() method. Returns a String object with its contents initialized to a subst...
OpenMS::PeptideIndexing::mm_max_
Int mm_max_
Definition: PeptideIndexing.h:843
OpenMS::ProteaseDigestion
Class for the enzymatic digestion of proteins.
Definition: ProteaseDigestion.h:60
SeqanIncludeWrapper.h
OpenMS::File::findDatabase
static String findDatabase(const String &db_name)
OpenMS::PeptideIndexing::DATABASE_EMPTY
Definition: PeptideIndexing.h:133
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::AABefore
char AABefore
Definition: PeptideIndexing.h:736
OpenMS::PeptideIndexing
Refreshes the protein references for all peptide hits in a vector of PeptideIdentifications and adds ...
Definition: PeptideIndexing.h:124
OpenMS::File::readable
static bool readable(const String &file)
Return true if the file exists and is readable.
OpenMS::PeptideEvidence::C_TERMINAL_AA
static const char C_TERMINAL_AA
Definition: PeptideEvidence.h:61
OpenMS::PeptideIndexing::allow_unmatched_
bool allow_unmatched_
Definition: PeptideIndexing.h:839
PeptideIndexing.h
OpenMS::PeptideIndexing::FoundProteinFunctor::xtandem_
bool xtandem_
Definition: PeptideIndexing.h:763
OpenMS::PeptideIndexing::enzyme_specificity_
String enzyme_specificity_
Definition: PeptideIndexing.h:834
OpenMS::Param
Management and storage of parameters / INI files.
Definition: Param.h:73
OpenMS::PeptideIndexing::FoundProteinFunctor
Definition: PeptideIndexing.h:753
OpenMS::String::remove
String & remove(char what)
Remove all occurrences of the character what.
seqan::isAmbiguous
bool isAmbiguous(AAcid c)
Definition: AhoCorasickAmbiguous.h:578
OpenMS::PeptideIndexing::IL_equivalent_
bool IL_equivalent_
Definition: PeptideIndexing.h:840
OpenMS::Map
Map class based on the STL map (containing several convenience functions)
Definition: Map.h:50
OPENMS_LOG_INFO
#define OPENMS_LOG_INFO
Macro if a information, e.g. a status should be reported.
Definition: LogStream.h:465
OpenMS::EnzymaticDigestion::setSpecificity
void setSpecificity(Specificity spec)
Sets the specificity for the digestion (default is SPEC_FULL).
OpenMS::StringUtils::toUpper
static String & toUpper(String &this_s)
Definition: StringUtils.h:832
OpenMS::PeptideIndexing::EXECUTION_OK
Definition: PeptideIndexing.h:132
PeptideIdentification.h
OpenMS::PeptideIndexing::ILLEGAL_PARAMETERS
Definition: PeptideIndexing.h:135
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::position
OpenMS::Int position
Definition: PeptideIndexing.h:735
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::tie
const std::tuple< const Size &, const Int &, const char &, const char & > tie() const
Definition: PeptideIndexing.h:739
seqan::position
Size< TNeedle >::Type position(const PatternAuxData< TNeedle > &dh)
Definition: AhoCorasickAmbiguous.h:561
StandardTypes.h
File.h
OpenMS::DecoyHelper::findDecoyString
static Result findDecoyString(FASTAContainer< T > &proteins)
Heuristic to determine the decoy string given a set of protein names.
Definition: FASTAContainer.h:359
OpenMS::ProgressLogger::setLogType
void setLogType(LogType type) const
Sets the progress log that should be used. The default type is NONE!
OpenMS::FASTAContainer
template parameter for vector-based FASTA access
Definition: FASTAContainer.h:82
TOPPBase.h
OpenMS::IdXMLFile
Used to load and store idXML files.
Definition: IdXMLFile.h:63
OpenMS::FASTAContainer< TFI_File >
FASTAContainer<TFI_File> will make FASTA entries available chunk-wise from start to end by loading it...
Definition: FASTAContainer.h:93
OpenMS::PeptideIndexing::keep_unreferenced_proteins_
bool keep_unreferenced_proteins_
Definition: PeptideIndexing.h:838