FASTAContainer.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Chris Bielow $
// $Authors: Chris Bielow $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/CONCEPT/LogStream.h>
#include <OpenMS/DATASTRUCTURES/String.h>
#include <OpenMS/DATASTRUCTURES/ListUtils.h>
#include <OpenMS/DATASTRUCTURES/StringUtilsSimple.h>
#include <OpenMS/FORMAT/FASTAFile.h>
 
#include <functional>
#include <fstream>
#include <unordered_map>
#include <memory>
#include <utility>
#include <vector>
 
#include <boost/regex.hpp>
 
namespace OpenMS
{
 
  struct TFI_File; 
  struct TFI_Vector; 
 
template<typename TBackend>
class FASTAContainer; // prototype
 
template<>
class FASTAContainer<TFI_File>
{
public:
  FASTAContainer() = delete;
 
  FASTAContainer(const String& FASTA_file)
    : f_(),
    offsets_(),
    data_fg_(),
    data_bg_(),
    chunk_offset_(0),
    filename_(FASTA_file)
  {
    f_.readStart(FASTA_file);
  }
 
  size_t getChunkOffset() const
  {
    return chunk_offset_;
  }
 
  bool activateCache()
  {
    chunk_offset_ += data_fg_.size();
    data_fg_.swap(data_bg_);
    data_bg_.clear(); // just in case someone calls activateCache() multiple times...
    return !data_fg_.empty();
  }
 
  bool cacheChunk(int suggested_size)
  {
    data_bg_.clear();
    data_bg_.reserve(suggested_size);
    FASTAFile::FASTAEntry p;
    for (int i = 0; i < suggested_size; ++i)
    {
      std::streampos spos = f_.position();
      if (!f_.readNext(p)) break;
      data_bg_.push_back(std::move(p));
      offsets_.push_back(spos);
    }
    return !data_bg_.empty();
  }
 
  size_t chunkSize() const
  {
    return data_fg_.size();
  }
 
  const FASTAFile::FASTAEntry& chunkAt(size_t pos) const
  {
    return data_fg_[pos];
  }
 
  bool readAt(FASTAFile::FASTAEntry& protein, size_t pos)
  {
    // check if position is currently cached...
    if (chunk_offset_ <= pos && pos < chunk_offset_ + chunkSize())
    {
      protein = data_fg_[pos - chunk_offset_];
      return true;
    }
    // read anew from disk...
    if (pos >= offsets_.size())
    {
      throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, pos, offsets_.size());
    }
    std::streampos spos = f_.position(); // save old position
    if (!f_.setPosition(offsets_[pos])) return false;
    bool r = f_.readNext(protein);
    f_.setPosition(spos); // restore old position
    return r;
  }
 
  bool empty()
  { // trusting the FASTA file can be read...
    return f_.atEnd() && offsets_.empty();
  }
 
  void reset()
  {
    offsets_.clear();
    data_fg_.clear();
    data_bg_.clear();
    chunk_offset_ = 0;
    f_.readStart(filename_);
  }
 
 
  size_t size() const
  {
    return offsets_.size();
  }
 
private:
  FASTAFile f_; 
  std::vector<std::streampos> offsets_; 
  std::vector<FASTAFile::FASTAEntry> data_fg_; 
  std::vector<FASTAFile::FASTAEntry> data_bg_; 
  size_t chunk_offset_; 
  std::string filename_;
};
 
template<>
class FASTAContainer<TFI_Vector>
{
public:
  FASTAContainer() = delete;
 
  FASTAContainer(const std::vector<FASTAFile::FASTAEntry>& data)
    : data_(data)
  {
  }
 
  size_t getChunkOffset() const
  {
    return 0;
  }
 
  bool activateCache()
  {
    if (!activate_count_)
    { 
      activate_count_ = 1;
      return true;
    }
    return false; 
  }
 
  bool cacheChunk(int /*suggested_size*/)
  {
    if (!cache_count_)
    { 
      cache_count_ = 1;
      return true;
    }
    return false; 
  }
 
  size_t chunkSize() const
  {
    return data_.size();
  }
 
  const FASTAFile::FASTAEntry& chunkAt(size_t pos) const
  {
    return data_[pos];
  }
 
  bool readAt(FASTAFile::FASTAEntry& protein, size_t pos) const
  {
    protein = data_[pos];
    return true;
  }
 
  bool empty() const
  {
    return data_.empty();
  }
 
  size_t size() const
  {
    return data_.size();
  }
 
  void reset()
  {
    activate_count_ = 0;
    cache_count_ = 0;
  }
 
private:
  const std::vector<FASTAFile::FASTAEntry>& data_; 
  int activate_count_ = 0;
  int cache_count_ = 0;
};
 
class DecoyHelper
{
public:
  struct Result
  {
    bool success; 
    String name; 
    bool is_prefix; 
 
    bool operator==(const Result& rhs) const
    {
      return success == rhs.success
             && name == rhs.name
             && is_prefix == rhs.is_prefix;
    }
  };
 
  struct DecoyStatistics
  {
    std::unordered_map<std::string, std::pair<Size, Size>> decoy_count; 
    std::unordered_map<std::string, std::string> decoy_case_sensitive; 
    Size all_prefix_occur{0}; 
    Size all_suffix_occur{0}; 
    Size all_proteins_count{0}; 
 
    bool operator==(const DecoyStatistics& rhs) const
    {
      return decoy_count == rhs.decoy_count
             && decoy_case_sensitive == rhs.decoy_case_sensitive
             && all_prefix_occur == rhs.all_prefix_occur
             && all_suffix_occur == rhs.all_suffix_occur
             && all_proteins_count == rhs.all_proteins_count;
    }
  };
 
  // decoy strings
  inline static const std::vector<std::string> affixes = { "decoy", "dec", "reverse", "rev", "reversed", "__id_decoy", "xxx", "shuffled", "shuffle", "pseudo", "random" };
 
  // setup prefix- and suffix regex strings
  inline static const std::string regexstr_prefix = std::string("^(") + ListUtils::concatenate<std::string>(affixes, "_*|") + "_*)";
  inline static const std::string regexstr_suffix = std::string("(_") + ListUtils::concatenate<std::string>(affixes, "*|_") + ")$";
 
  template<typename T>
  static Result findDecoyString(FASTAContainer<T>& proteins)
  {
    // calls function to search for decoys in input data
    DecoyStatistics decoy_stats = countDecoys(proteins);
 
    // DEBUG ONLY: print counts of found decoys
    for (const auto &a : decoy_stats.decoy_count)
    {
      OPENMS_LOG_DEBUG << a.first << "\t" << a.second.first << "\t" << a.second.second << std::endl;
    }
 
    // less than 40% of proteins are decoys -> won't be able to determine a decoy string and its position
    // return default values
    if (static_cast<double>(decoy_stats.all_prefix_occur + decoy_stats.all_suffix_occur) < 0.4 * static_cast<double>(decoy_stats.all_proteins_count))
    {
      OPENMS_LOG_ERROR << "Unable to determine decoy string (not enough occurrences; <40%)!" << std::endl;
      return {false, "?", true};
    }
 
    if (decoy_stats.all_prefix_occur == decoy_stats.all_suffix_occur)
    {
      OPENMS_LOG_ERROR << "Unable to determine decoy string (prefix and suffix occur equally often)!" << std::endl;
      return {false, "?", true};
    }
 
    // Decoy prefix occurred at least 80% of all prefixes + observed in at least 40% of all proteins -> set it as prefix decoy
    for (const auto& pair : decoy_stats.decoy_count)
    {
      const std::string & case_insensitive_decoy_string = pair.first;
      const std::pair<Size, Size>& prefix_suffix_counts = pair.second;
      double freq_prefix = static_cast<double>(prefix_suffix_counts.first) / static_cast<double>(decoy_stats.all_prefix_occur);
      double freq_prefix_in_proteins = static_cast<double>(prefix_suffix_counts.first) / static_cast<double>(decoy_stats.all_proteins_count);
 
      if (freq_prefix >= 0.8 && freq_prefix_in_proteins >= 0.4)
      {
        if (prefix_suffix_counts.first != decoy_stats.all_prefix_occur)
        {
          OPENMS_LOG_WARN << "More than one decoy prefix observed!" << std::endl;
          OPENMS_LOG_WARN << "Using most frequent decoy prefix (" << (int)(freq_prefix * 100) << "%)" << std::endl;
        }
 
        return { true, decoy_stats.decoy_case_sensitive[case_insensitive_decoy_string], true};
      }
    }
 
    // Decoy suffix occurred at least 80% of all suffixes + observed in at least 40% of all proteins -> set it as suffix decoy
    for (const auto& pair : decoy_stats.decoy_count)
    {
      const std::string& case_insensitive_decoy_string = pair.first;
      const std::pair<Size, Size>& prefix_suffix_counts = pair.second;
      double freq_suffix = static_cast<double>(prefix_suffix_counts.second) / static_cast<double>(decoy_stats.all_suffix_occur);
      double freq_suffix_in_proteins = static_cast<double>(prefix_suffix_counts.second) / static_cast<double>(decoy_stats.all_proteins_count);
 
      if (freq_suffix >= 0.8 && freq_suffix_in_proteins >= 0.4)
      {
        if (prefix_suffix_counts.second != decoy_stats.all_suffix_occur)
        {
          OPENMS_LOG_WARN << "More than one decoy suffix observed!" << std::endl;
          OPENMS_LOG_WARN << "Using most frequent decoy suffix (" << (int)(freq_suffix * 100) << "%)" << std::endl;
        }
 
        return { true, decoy_stats.decoy_case_sensitive[case_insensitive_decoy_string], false};
      }
    }
 
    OPENMS_LOG_ERROR << "Unable to determine decoy string and its position. Please provide a decoy string and its position as parameters." << std::endl;
    return {false, "?", true};
  }
 
  template<typename T>
  static DecoyStatistics countDecoys(FASTAContainer<T>& proteins)
  {
    // common decoy strings in FASTA files
    // note: decoy prefixes/suffices must be provided in lower case
 
    DecoyStatistics ds;
 
    // setup regexes
    const boost::regex pattern_prefix(regexstr_prefix);
    const boost::regex pattern_suffix(regexstr_suffix);
 
    constexpr size_t PROTEIN_CACHE_SIZE = 4e5;
 
    while (true)
    {
      proteins.cacheChunk(PROTEIN_CACHE_SIZE);
      if (!proteins.activateCache()) break;
 
      auto prot_count = (SignedSize)proteins.chunkSize();
      ds.all_proteins_count += prot_count;
 
      boost::smatch sm;
      for (SignedSize i = 0; i < prot_count; ++i)
      {
        String seq = proteins.chunkAt(i).identifier;
 
        String seq_lower = seq;
        seq_lower.toLower();
 
        // search for prefix
        bool found_prefix = boost::regex_search(seq_lower, sm, pattern_prefix);
        if (found_prefix)
        {
          std::string match = sm[0];
          ds.all_prefix_occur++;
 
          // increase count of observed prefix
          ds.decoy_count[match].first++;
 
          // store observed (case sensitive and with special characters)
          std::string seq_decoy = StringUtils::prefix(seq, match.length());
          ds.decoy_case_sensitive[match] = seq_decoy;
        }
 
        // search for suffix
        bool found_suffix = boost::regex_search(seq_lower, sm, pattern_suffix);
        if (found_suffix)
        {
          std::string match = sm[0];
          ds.all_suffix_occur++;
 
          // increase count of observed suffix
          ds.decoy_count[match].second++;
 
          // store observed (case sensitive and with special characters)
          std::string seq_decoy = StringUtils::suffix(seq, match.length());
          ds.decoy_case_sensitive[match] = seq_decoy;
        }
      }
    }
    return ds;
  }
 
private:
  using DecoyStringToAffixCount = std::unordered_map<std::string, std::pair<Size, Size>>;
  using CaseInsensitiveToCaseSensitiveDecoy = std::unordered_map<std::string, std::string>;
};
 
} // namespace OpenMS