OpenMS  2.8.0
FASTAContainer.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-2021.
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: Chris Bielow $
33 // --------------------------------------------------------------------------
34 
35 #pragma once
36 
43 
44 #include <functional>
45 #include <fstream>
46 #include <unordered_map>
47 #include <memory>
48 #include <utility>
49 #include <vector>
50 
51 #include <boost/regex.hpp>
52 
53 namespace OpenMS
54 {
55 
56  struct TFI_File;
57  struct TFI_Vector;
58 
81 template<typename TBackend>
82 class FASTAContainer; // prototype
83 
92 template<>
93 class FASTAContainer<TFI_File>
94 {
95 public:
96  FASTAContainer() = delete;
97 
99  FASTAContainer(const String& FASTA_file)
100  : f_(),
101  offsets_(),
102  data_fg_(),
103  data_bg_(),
104  chunk_offset_(0),
105  filename_(FASTA_file)
106  {
107  f_.readStart(FASTA_file);
108  }
109 
111  size_t getChunkOffset() const
112  {
113  return chunk_offset_;
114  }
115 
123  {
124  chunk_offset_ += data_fg_.size();
125  data_fg_.swap(data_bg_);
126  data_bg_.clear(); // just in case someone calls activateCache() multiple times...
127  return !data_fg_.empty();
128  }
129 
136  bool cacheChunk(int suggested_size)
137  {
138  data_bg_.clear();
139  data_bg_.reserve(suggested_size);
141  for (int i = 0; i < suggested_size; ++i)
142  {
143  std::streampos spos = f_.position();
144  if (!f_.readNext(p)) break;
145  data_bg_.push_back(std::move(p));
146  offsets_.push_back(spos);
147  }
148  return !data_bg_.empty();
149  }
150 
152  size_t chunkSize() const
153  {
154  return data_fg_.size();
155  }
156 
164  const FASTAFile::FASTAEntry& chunkAt(size_t pos) const
165  {
166  return data_fg_[pos];
167  }
168 
182  bool readAt(FASTAFile::FASTAEntry& protein, size_t pos)
183  {
184  // check if position is currently cached...
185  if (chunk_offset_ <= pos && pos < chunk_offset_ + chunkSize())
186  {
187  protein = data_fg_[pos - chunk_offset_];
188  return true;
189  }
190  // read anew from disk...
191  if (pos >= offsets_.size())
192  {
193  throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, pos, offsets_.size());
194  }
195  std::streampos spos = f_.position(); // save old position
196  if (!f_.setPosition(offsets_[pos])) return false;
197  bool r = f_.readNext(protein);
198  f_.setPosition(spos); // restore old position
199  return r;
200  }
201 
203  bool empty()
204  { // trusting the FASTA file can be read...
205  return f_.atEnd() && offsets_.empty();
206  }
207 
209  void reset()
210  {
211  offsets_.clear();
212  data_fg_.clear();
213  data_bg_.clear();
214  chunk_offset_ = 0;
215  f_.readStart(filename_);
216  }
217 
218 
224  size_t size() const
225  {
226  return offsets_.size();
227  }
228 
229 private:
231  std::vector<std::streampos> offsets_;
232  std::vector<FASTAFile::FASTAEntry> data_fg_;
233  std::vector<FASTAFile::FASTAEntry> data_bg_;
234  size_t chunk_offset_;
235  std::string filename_;
236 };
237 
244 template<>
245 class FASTAContainer<TFI_Vector>
246 {
247 public:
248  FASTAContainer() = delete;
249 
255  FASTAContainer(const std::vector<FASTAFile::FASTAEntry>& data)
256  : data_(data)
257  {
258  }
259 
261  size_t getChunkOffset() const
262  {
263  return 0;
264  }
265 
271  {
272  if (!activate_count_)
273  {
274  activate_count_ = 1;
275  return true;
276  }
277  return false;
278  }
279 
283  bool cacheChunk(int /*suggested_size*/)
284  {
285  if (!cache_count_)
286  {
287  cache_count_ = 1;
288  return true;
289  }
290  return false;
291  }
292 
297  size_t chunkSize() const
298  {
299  return data_.size();
300  }
301 
303  const FASTAFile::FASTAEntry& chunkAt(size_t pos) const
304  {
305  return data_[pos];
306  }
307 
309  bool readAt(FASTAFile::FASTAEntry& protein, size_t pos) const
310  {
311  protein = data_[pos];
312  return true;
313  }
314 
316  bool empty() const
317  {
318  return data_.empty();
319  }
320 
322  size_t size() const
323  {
324  return data_.size();
325  }
326 
328  void reset()
329  {
330  activate_count_ = 0;
331  cache_count_ = 0;
332  }
333 
334 private:
335  const std::vector<FASTAFile::FASTAEntry>& data_;
336  int activate_count_ = 0;
337  int cache_count_ = 0;
338 };
339 
344 {
345 public:
346  struct Result
347  {
348  bool success;
350  bool is_prefix;
351  };
352 
353  // decoy strings
354  inline static const std::vector<std::string> affixes = { "decoy", "dec", "reverse", "rev", "reversed", "__id_decoy", "xxx", "shuffled", "shuffle", "pseudo", "random" };
355 
356  // setup prefix- and suffix regex strings
357  inline static const std::string regexstr_prefix = std::string("^(") + ListUtils::concatenate<std::string>(affixes, "_*|") + "_*)";
358  inline static const std::string regexstr_suffix = std::string("(_") + ListUtils::concatenate<std::string>(affixes, "*|_") + ")$";
359 
367  template<typename T>
369  {
370  // common decoy strings in FASTA files
371  // note: decoy prefixes/suffices must be provided in lower case
372 
373  // map decoys to counts of occurrences as prefix/suffix
374  DecoyStringToAffixCount decoy_count;
375  // map case insensitive strings back to original case (as used in fasta)
376  CaseInsensitiveToCaseSensitiveDecoy decoy_case_sensitive;
377 
378  // setup regexes
379  const boost::regex pattern_prefix(regexstr_prefix);
380  const boost::regex pattern_suffix(regexstr_suffix);
381 
382  Size all_prefix_occur(0), all_suffix_occur(0), all_proteins_count(0);
383 
384  constexpr size_t PROTEIN_CACHE_SIZE = 4e5;
385 
386  while (true)
387  {
388  proteins.cacheChunk(PROTEIN_CACHE_SIZE);
389  if (!proteins.activateCache()) break;
390 
391  auto prot_count = (SignedSize)proteins.chunkSize();
392  all_proteins_count += prot_count;
393 
394  boost::smatch sm;
395  for (SignedSize i = 0; i < prot_count; ++i)
396  {
397  String seq = proteins.chunkAt(i).identifier;
398 
399  String seq_lower = seq;
400  seq_lower.toLower();
401 
402  // search for prefix
403  bool found_prefix = boost::regex_search(seq_lower, sm, pattern_prefix);
404  if (found_prefix)
405  {
406  std::string match = sm[0];
407  all_prefix_occur++;
408 
409  // increase count of observed prefix
410  decoy_count[match].first++;
411 
412  // store observed (case sensitive and with special characters)
413  std::string seq_decoy = StringUtils::prefix(seq, match.length());
414  decoy_case_sensitive[match] = seq_decoy;
415  }
416 
417  // search for suffix
418  bool found_suffix = boost::regex_search(seq_lower, sm, pattern_suffix);
419  if (found_suffix)
420  {
421  std::string match = sm[0];
422  all_suffix_occur++;
423 
424  // increase count of observed suffix
425  decoy_count[match].second++;
426 
427  // store observed (case sensitive and with special characters)
428  std::string seq_decoy = StringUtils::suffix(seq, match.length());
429  decoy_case_sensitive[match] = seq_decoy;
430  }
431  }
432  }
433 
434  // DEBUG ONLY: print counts of found decoys
435  for (auto &a : decoy_count)
436  {
437  OPENMS_LOG_DEBUG << a.first << "\t" << a.second.first << "\t" << a.second.second << std::endl;
438  }
439 
440  // less than 40% of proteins are decoys -> won't be able to determine a decoy string and its position
441  // return default values
442  if (static_cast<double>(all_prefix_occur + all_suffix_occur) < 0.4 * static_cast<double>(all_proteins_count))
443  {
444  OPENMS_LOG_ERROR << "Unable to determine decoy string (not enough occurrences; <40%)!" << std::endl;
445  return {false, "?", true};
446  }
447 
448  if (all_prefix_occur == all_suffix_occur)
449  {
450  OPENMS_LOG_ERROR << "Unable to determine decoy string (prefix and suffix occur equally often)!" << std::endl;
451  return {false, "?", true};
452  }
453 
454  // Decoy prefix occurred at least 80% of all prefixes + observed in at least 40% of all proteins -> set it as prefix decoy
455  for (const auto& pair : decoy_count)
456  {
457  const std::string & case_insensitive_decoy_string = pair.first;
458  const std::pair<Size, Size>& prefix_suffix_counts = pair.second;
459  double freq_prefix = static_cast<double>(prefix_suffix_counts.first) / static_cast<double>(all_prefix_occur);
460  double freq_prefix_in_proteins = static_cast<double>(prefix_suffix_counts.first) / static_cast<double>(all_proteins_count);
461 
462  if (freq_prefix >= 0.8 && freq_prefix_in_proteins >= 0.4)
463  {
464  if (prefix_suffix_counts.first != all_prefix_occur)
465  {
466  OPENMS_LOG_WARN << "More than one decoy prefix observed!" << std::endl;
467  OPENMS_LOG_WARN << "Using most frequent decoy prefix (" << (int)(freq_prefix * 100) << "%)" << std::endl;
468  }
469 
470  return { true, decoy_case_sensitive[case_insensitive_decoy_string], true};
471  }
472  }
473 
474  // Decoy suffix occurred at least 80% of all suffixes + observed in at least 40% of all proteins -> set it as suffix decoy
475  for (const auto& pair : decoy_count)
476  {
477  const std::string& case_insensitive_decoy_string = pair.first;
478  const std::pair<Size, Size>& prefix_suffix_counts = pair.second;
479  double freq_suffix = static_cast<double>(prefix_suffix_counts.second) / static_cast<double>(all_suffix_occur);
480  double freq_suffix_in_proteins = static_cast<double>(prefix_suffix_counts.second) / static_cast<double>(all_proteins_count);
481 
482  if (freq_suffix >= 0.8 && freq_suffix_in_proteins >= 0.4)
483  {
484  if (prefix_suffix_counts.second != all_suffix_occur)
485  {
486  OPENMS_LOG_WARN << "More than one decoy suffix observed!" << std::endl;
487  OPENMS_LOG_WARN << "Using most frequent decoy suffix (" << (int)(freq_suffix * 100) << "%)" << std::endl;
488  }
489 
490  return { true, decoy_case_sensitive[case_insensitive_decoy_string], false};
491  }
492  }
493 
494  OPENMS_LOG_ERROR << "Unable to determine decoy string and its position. Please provide a decoy string and its position as parameters." << std::endl;
495  return {false, "?", true};
496  }
497 
498 private:
499  using DecoyStringToAffixCount = std::unordered_map<std::string, std::pair<Size, Size>>;
500  using CaseInsensitiveToCaseSensitiveDecoy = std::unordered_map<std::string, std::string>;
501 };
502 
503 } // namespace OpenMS
504 
#define OPENMS_LOG_DEBUG
Macro for general debugging information.
Definition: LogStream.h:470
#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
#define OPENMS_LOG_ERROR
Macro to be used if non-fatal error are reported (processing continues)
Definition: LogStream.h:455
Helper class for calculations on decoy proteins.
Definition: FASTAContainer.h:344
std::unordered_map< std::string, std::pair< Size, Size > > DecoyStringToAffixCount
Definition: FASTAContainer.h:499
std::unordered_map< std::string, std::string > CaseInsensitiveToCaseSensitiveDecoy
Definition: FASTAContainer.h:500
bool is_prefix
on success, was it a prefix or suffix
Definition: FASTAContainer.h:350
static Result findDecoyString(FASTAContainer< T > &proteins)
Heuristic to determine the decoy string given a set of protein names.
Definition: FASTAContainer.h:368
static const std::string regexstr_prefix
Definition: FASTAContainer.h:357
bool success
did more than 40% of proteins have the *same* prefix or suffix
Definition: FASTAContainer.h:348
String name
on success, what was the decoy string?
Definition: FASTAContainer.h:349
static const std::vector< std::string > affixes
Definition: FASTAContainer.h:354
static const std::string regexstr_suffix
Definition: FASTAContainer.h:358
Definition: FASTAContainer.h:347
Int overflow exception.
Definition: Exception.h:247
std::vector< FASTAFile::FASTAEntry > data_bg_
prefetched (background) data; will become the next active data
Definition: FASTAContainer.h:233
bool readAt(FASTAFile::FASTAEntry &protein, size_t pos)
Retrieve a FASTA entry at global position pos (must not be behind the currently active chunk,...
Definition: FASTAContainer.h:182
std::string filename_
FASTA file name.
Definition: FASTAContainer.h:235
size_t size() const
NOT the number of entries in the FASTA file, but merely the number of already read entries (since we ...
Definition: FASTAContainer.h:224
bool empty()
is the FASTA file empty?
Definition: FASTAContainer.h:203
std::vector< std::streampos > offsets_
internal byte offsets into FASTA file for random access reading of previous entries.
Definition: FASTAContainer.h:231
bool activateCache()
Swaps in the background cache of entries, read previously via cacheChunk()
Definition: FASTAContainer.h:122
FASTAContainer(const String &FASTA_file)
C'tor with FASTA filename.
Definition: FASTAContainer.h:99
size_t chunk_offset_
number of entries before the current chunk
Definition: FASTAContainer.h:234
bool cacheChunk(int suggested_size)
Prefetch a new cache in the background, with up to suggested_size entries (or fewer upon reaching end...
Definition: FASTAContainer.h:136
FASTAFile f_
FASTA file connection.
Definition: FASTAContainer.h:230
std::vector< FASTAFile::FASTAEntry > data_fg_
active (foreground) data
Definition: FASTAContainer.h:232
size_t chunkSize() const
number of entries in active cache
Definition: FASTAContainer.h:152
void reset()
resets reading of the FASTA file, enables fresh reading of the FASTA from the beginning
Definition: FASTAContainer.h:209
const FASTAFile::FASTAEntry & chunkAt(size_t pos) const
Retrieve a FASTA entry at cache position pos (fast)
Definition: FASTAContainer.h:164
size_t getChunkOffset() const
how many entries were read and got swapped out already
Definition: FASTAContainer.h:111
FASTAContainer(const std::vector< FASTAFile::FASTAEntry > &data)
C'tor for already existing data (by reference).
Definition: FASTAContainer.h:255
size_t size() const
calls size() on underlying vector
Definition: FASTAContainer.h:322
bool readAt(FASTAFile::FASTAEntry &protein, size_t pos) const
fast access to an entry
Definition: FASTAContainer.h:309
bool empty() const
calls empty() on underlying vector
Definition: FASTAContainer.h:316
bool cacheChunk(int)
no-op (since data is already fully available as vector)
Definition: FASTAContainer.h:283
bool activateCache()
no-op (since data is already fully available as vector)
Definition: FASTAContainer.h:270
const std::vector< FASTAFile::FASTAEntry > & data_
reference to existing data
Definition: FASTAContainer.h:335
size_t chunkSize() const
active data spans the full range, i.e. size of container
Definition: FASTAContainer.h:297
void reset()
required for template parameters!
Definition: FASTAContainer.h:328
const FASTAFile::FASTAEntry & chunkAt(size_t pos) const
fast access to chunked (i.e. all) entries
Definition: FASTAContainer.h:303
size_t getChunkOffset() const
always 0, since this specialization requires/supports no chunking
Definition: FASTAContainer.h:261
This class serves for reading in and writing FASTA files If the protein/gene sequence contains unusua...
Definition: FASTAFile.h:61
A more convenient string class.
Definition: String.h:60
String & toLower()
Converts the string to lowercase.
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
static String suffix(const String &this_s, size_t length)
Definition: StringUtilsSimple.h:156
static String prefix(const String &this_s, size_t length)
Definition: StringUtilsSimple.h:147
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
template parameter for vector-based FASTA access
Definition: FASTAContainer.h:82
FASTA entry type (identifier, description and sequence) The first String corresponds to the identifie...
Definition: FASTAFile.h:72