| 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
38 #include <OpenMS/ANALYSIS/ID/AhoCorasickAmbiguous.h>
39 #include <OpenMS/CHEMISTRY/ProteaseDigestion.h>
40 #include <OpenMS/CHEMISTRY/ProteaseDB.h>
41 #include <OpenMS/CONCEPT/LogStream.h>
42 #include <OpenMS/CONCEPT/ProgressLogger.h>
43 #include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
44 #include <OpenMS/DATASTRUCTURES/FASTAContainer.h>
45 #include <OpenMS/DATASTRUCTURES/ListUtils.h>
46 #include <OpenMS/DATASTRUCTURES/StringUtils.h>
47 #include <OpenMS/DATASTRUCTURES/SeqanIncludeWrapper.h>
48 #include <OpenMS/FORMAT/FASTAFile.h>
49 #include <OpenMS/KERNEL/StandardTypes.h>
50 #include <OpenMS/METADATA/PeptideEvidence.h>
51 #include <OpenMS/METADATA/PeptideIdentification.h>
52 #include <OpenMS/METADATA/ProteinIdentification.h>
53 #include <OpenMS/SYSTEM/StopWatch.h>
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
130 enum ExitCodes
131 {
132 EXECUTION_OK,
133 DATABASE_EMPTY,
134 PEPTIDE_IDS_EMPTY,
135 ILLEGAL_PARAMETERS,
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);
279 AhoCorasickAmbiguous::PeptideDB pep_DB;
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;
321 SysInfo::MemUsage mu;
322 OPENMS_LOG_INFO << "Building trie ...";
323 StopWatch s;
324 s.start();
325 AhoCorasickAmbiguous::FuzzyACPattern pattern;
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
627 FASTAFile::FASTAEntry fe;
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
732 struct PeptideProteinMatchInformation
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 }
747 bool operator==(const PeptideProteinMatchInformation& other) const
748 {
749 return tie() == other.tie();
750 }
751 };
752
753 struct FoundProteinFunctor
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:
762 ProteaseDigestion enzyme_;
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
771 void merge(FoundProteinFunctor& other)
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,
796 OpenMS::Int position)
797 {
798 if (enzyme_.isValidProduct(seq_prot, position, len_pep, true, true, xtandem_))
799 {
800 PeptideProteinMatchInformation match;
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
830 String decoy_string_;
831 bool prefix_;
832 String missing_decoy_action_;
833 String enzyme_name_;
834 String enzyme_specificity_;
835
836 bool write_protein_sequence_;
837 bool write_protein_description_;
838 bool keep_unreferenced_proteins_;
839 bool allow_unmatched_;
840 bool IL_equivalent_;
841
842 Int aaa_max_;
843 Int mm_max_;
844 };
845 }
846
OpenMS::StopWatch::getClockTime
double getClockTime() const
OpenMS::SignedSize
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
OpenMS::SysInfo::MemUsage
A convenience class to report either absolute or delta (between two timepoints) RAM usage...
Definition: SysInfo.h:83
ProteinIdentification.h
OpenMS::PeptideIndexing::prefix_
bool prefix_
Definition: PeptideIndexing.h:831
OpenMS::String
A more convenient string class.
Definition: String.h:58
OpenMS::PeptideIndexing::DATABASE_EMPTY
Definition: PeptideIndexing.h:133
OpenMS::PeptideIndexing::PEPTIDE_IDS_EMPTY
Definition: PeptideIndexing.h:134
OpenMS::PeptideIndexing::EXECUTION_OK
Definition: PeptideIndexing.h:132
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::tie
const std::tuple< const Size &, const Int &, const char &, const char & > tie() const
Definition: PeptideIndexing.h:739
OpenMS::EnzymaticDigestion::setSpecificity
void setSpecificity(Specificity spec)
Sets the specificity for the digestion (default is SPEC_FULL).
OpenMS::String::hasSuffix
bool hasSuffix(const String &string) const
true if String ends with string, false otherwise
seqan::AAString
String< AAcid, Alloc< void > > AAString
Definition: AhoCorasickAmbiguous.h:206
OpenMS::PeptideIndexing::ExitCodes
ExitCodes
Exit codes.
Definition: PeptideIndexing.h:130
OpenMS::FASTAFile::FASTAEntry::sequence
String sequence
Definition: FASTAFile.h:80
OpenMS::DefaultParamHandler
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:91
DEBUG_ONLY
#define DEBUG_ONLY
Definition: AhoCorasickAmbiguous.h:46
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::operator<
bool operator<(const PeptideProteinMatchInformation &other) const
Definition: PeptideIndexing.h:743
ProteaseDigestion.h
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
StringUtils.h
seqan::position
Size< TNeedle >::Type position(const PatternAuxData< TNeedle > &dh)
Definition: AhoCorasickAmbiguous.h:561
OpenMS::PeptideIndexing::mm_max_
Int mm_max_
Definition: PeptideIndexing.h:843
OpenMS::String::substitute
String & substitute(char from, char to)
Replaces all occurrences of the character from by the character to.
OpenMS::AhoCorasickAmbiguous::getHitDBIndex
Size getHitDBIndex()
Get index of hit into peptide database of the pattern.
Definition: AhoCorasickAmbiguous.h:1047
OpenMS::Map::has
bool has(const Key &key) const
Test whether the map contains the given key.
Definition: Map.h:108
OpenMS::String::hasPrefix
bool hasPrefix(const String &string) const
true if String begins with string, false otherwise
OpenMS::PeptideIndexing::FoundProteinFunctor::filter_passed
OpenMS::Size filter_passed
Definition: PeptideIndexing.h:758
OpenMS::ProteinHit
Representation of a protein hit.
Definition: ProteinHit.h:57
StandardTypes.h
SeqanIncludeWrapper.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::PeptideEvidence
Representation of a peptide evidence.
Definition: PeptideEvidence.h:50
OpenMS::AhoCorasickAmbiguous::FuzzyACPattern
::seqan::Pattern< PeptideDB, ::seqan::FuzzyAC > FuzzyACPattern
Definition: AhoCorasickAmbiguous.h:974
OpenMS::PeptideIndexing::keep_unreferenced_proteins_
bool keep_unreferenced_proteins_
Definition: PeptideIndexing.h:838
OpenMS::PeptideEvidence::C_TERMINAL_AA
static const char C_TERMINAL_AA
Definition: PeptideEvidence.h:61
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::enzyme_name_
String enzyme_name_
Definition: PeptideIndexing.h:833
OpenMS::ProteaseDigestion
Class for the enzymatic digestion of proteins.
Definition: ProteaseDigestion.h:60
OpenMS::ProteinHit::setSequence
void setSequence(const String &sequence)
sets the protein sequence
OpenMS::ProteinHit::setDescription
void setDescription(const String &description)
sets the description of the protein
FASTAContainer.h
OpenMS::PeptideIndexing::decoy_string_
String decoy_string_
Definition: PeptideIndexing.h:830
OpenMS::PeptideIndexing
Refreshes the protein references for all peptide hits in a vector of PeptideIdentifications and adds ...
Definition: PeptideIndexing.h:124
AhoCorasickAmbiguous.h
OpenMS::String::has
bool has(Byte byte) const
true if String contains the byte, false otherwise
OpenMS::PeptideIndexing::FoundProteinFunctor::merge
void merge(FoundProteinFunctor &other)
Definition: PeptideIndexing.h:771
OpenMS::PeptideIndexing::IL_equivalent_
bool IL_equivalent_
Definition: PeptideIndexing.h:840
PeptideEvidence.h
ProteaseDB.h
OpenMS::FASTAContainer
template parameter for vector-based FASTA access
Definition: FASTAContainer.h:82
OpenMS::StopWatch::stop
bool stop()
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::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::PeptideIndexing::PeptideProteinMatchInformation::operator==
bool operator==(const PeptideProteinMatchInformation &other) const
Definition: PeptideIndexing.h:747
OpenMS::StopWatch
StopWatch Class.
Definition: StopWatch.h:59
OpenMS::AhoCorasickAmbiguous::findNext
bool findNext(const FuzzyACPattern &pattern)
Enumerate hits.
Definition: AhoCorasickAmbiguous.h:1037
OpenMS::String::remove
String & remove(char what)
Remove all occurrences of the character what.
FASTAFile.h
OpenMS::PeptideIndexing::FoundProteinFunctor
Definition: PeptideIndexing.h:753
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...
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
ProgressLogger.h
OpenMS::PeptideIndexing::FoundProteinFunctor::xtandem_
bool xtandem_
Definition: PeptideIndexing.h:763
OpenMS::Map
Map class based on the STL map (containing several convenience functions)
Definition: Map.h:50
OpenMS::FASTAContainer< TFI_Vector >
FASTAContainer<TFI_Vector> simply takes an existing vector of FASTAEntries and provides the same inte...
Definition: FASTAContainer.h:243
OpenMS::FASTAFile::FASTAEntry::description
String description
Definition: FASTAFile.h:79
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::position
OpenMS::Int position
Definition: PeptideIndexing.h:735
OpenMS::PeptideIndexing::missing_decoy_action_
String missing_decoy_action_
Definition: PeptideIndexing.h:832
OpenMS::PeptideIndexing::FoundProteinFunctor::FoundProteinFunctor
FoundProteinFunctor(const ProteaseDigestion &enzyme, bool xtandem)
Definition: PeptideIndexing.h:766
PeptideIdentification.h
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
OpenMS::PeptideIndexing::write_protein_sequence_
bool write_protein_sequence_
Definition: PeptideIndexing.h:836
OpenMS::AhoCorasickAmbiguous::getHitProteinPosition
Int getHitProteinPosition()
Offset into protein sequence where hit was found.
Definition: AhoCorasickAmbiguous.h:1057
OpenMS::ProgressLogger
Base class for all classes that want to report their progress.
Definition: ProgressLogger.h:54
OPENMS_LOG_INFO
#define OPENMS_LOG_INFO
Macro if a information, e.g. a status should be reported.
Definition: LogStream.h:465
OpenMS::StopWatch::reset
void reset()
LogStream.h
OpenMS::PeptideIndexing::write_protein_description_
bool write_protein_description_
Definition: PeptideIndexing.h:837
OpenMS::FASTAFile::FASTAEntry
FASTA entry type (identifier, description and sequence)
Definition: FASTAFile.h:76
OpenMS::AhoCorasickAmbiguous
Extended Aho-Corasick algorithm capable of matching ambiguous amino acids in the pattern (i...
Definition: AhoCorasickAmbiguous.h:970
DefaultParamHandler.h
OpenMS::PeptideIndexing::PeptideProteinMatchInformation
Definition: PeptideIndexing.h:732
OpenMS::StopWatch::start
bool start()
OpenMS::SysInfo::MemUsage::delta
String delta(const String &event="delta")
OpenMS::EnzymaticDigestion::getSpecificityByName
static Specificity getSpecificityByName(const String &name)
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::EnzymaticDigestion::getEnzymeName
String getEnzymeName() const
Returns the enzyme for the digestion.
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::protein_index
OpenMS::Size protein_index
Definition: PeptideIndexing.h:734
OpenMS::PeptideEvidence::N_TERMINAL_AA
static const char N_TERMINAL_AA
Definition: PeptideEvidence.h:60
OpenMS::PeptideIndexing::aaa_max_
Int aaa_max_
Definition: PeptideIndexing.h:842
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::AABefore
char AABefore
Definition: PeptideIndexing.h:736
OpenMS::PeptideIndexing::allow_unmatched_
bool allow_unmatched_
Definition: PeptideIndexing.h:839
int
OpenMS::PeptideIndexing::enzyme_specificity_
String enzyme_specificity_
Definition: PeptideIndexing.h:834
OpenMS::AhoCorasickAmbiguous::PeptideDB
::seqan::StringSet<::seqan::AAString > PeptideDB
Definition: AhoCorasickAmbiguous.h:973
seqan::isAmbiguous
bool isAmbiguous(AAcid c)
Definition: AhoCorasickAmbiguous.h:578
OpenMS::PeptideIndexing::PeptideProteinMatchInformation::AAAfter
char AAAfter
Definition: PeptideIndexing.h:737
OpenMS::PeptideIndexing::FoundProteinFunctor::filter_rejected
OpenMS::Size filter_rejected
Definition: PeptideIndexing.h:759
OpenMS::PeptideIndexing::FoundProteinFunctor::pep_to_prot
MapType pep_to_prot
Definition: PeptideIndexing.h:757
OpenMS::MetaInfoInterface::setMetaValue
void setMetaValue(const String &name, const DataValue &value)
Sets the DataValue corresponding to a name.
OPENMS_LOG_ERROR
#define OPENMS_LOG_ERROR
Macro to be used if non-fatal error are reported (processing continues)
Definition: LogStream.h:455
ListUtils.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::ProteinHit::setAccession
void setAccession(const String &accession)
sets the accession of the protein
OpenMS::ProteaseDigestion::setEnzyme
void setEnzyme(const String &name)
Sets the enzyme for the digestion (by name)
OpenMS::StringUtils::toUpper
static String & toUpper(String &this_s)
Definition: StringUtils.h:832
OpenMS::PeptideIndexing::FoundProteinFunctor::MapType
std::map< OpenMS::Size, std::set< PeptideProteinMatchInformation > > MapType
Definition: PeptideIndexing.h:756
SysInfo.h
OpenMS::StopWatch::toString
String toString() const
get a compact representation of the current time status.
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::SysInfo::MemUsage::after
void after()
record data for the second timepoint
OpenMS::PeptideIndexing::ILLEGAL_PARAMETERS
Definition: PeptideIndexing.h:135
OpenMS::Size
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
OpenMS::PeptideIndexing::FoundProteinFunctor::enzyme_
ProteaseDigestion enzyme_
Definition: PeptideIndexing.h:762
StopWatch.h