FASTAContainer.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2020.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $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/StringUtils.h>
#include <OpenMS/FORMAT/FASTAFile.h>
 
#include <functional>
#include <fstream>
#include <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)
  {
    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() const
  { // trusting the FASTA file can be read...
    return f_.atEnd() && offsets_.empty();
  }
 
  void reset()
  {
    f_.setPosition(0);
    offsets_.clear();
    data_fg_.clear();
    data_bg_.clear();
    chunk_offset_ = 0;
  }
 
 
  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_; 
};
 
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; //< did >=40% of proteins have the *same* prefix or suffix
    String name; //< on success, what was the decoy string?
    bool is_prefix; //< on success, was it a prefix or suffix
  };
 
  template<typename T>
  static Result findDecoyString(FASTAContainer<T>& proteins)
  {
    // common decoy strings in FASTA files
    // note: decoy prefixes/suffices must be provided in lower case
    const std::vector<std::string> affixes{ "decoy", "dec", "reverse", "rev", "reversed", "__id_decoy", "xxx", "shuffled", "shuffle", "pseudo", "random" };
 
    // map decoys to counts of occurrences as prefix/suffix
    DecoyStringToAffixCount decoy_count;
    // map case insensitive strings back to original case (as used in fasta)
    CaseInsensitiveToCaseSensitiveDecoy decoy_case_sensitive;
 
    // setup prefix- and suffix regex strings
    const std::string regexstr_prefix = std::string("^(") + ListUtils::concatenate<std::string>(affixes, "_*|") + "_*)";
    const std::string regexstr_suffix = std::string("(") + ListUtils::concatenate<std::string>(affixes, "_*|") + "_*)$";
 
    // setup regexes
    const boost::regex pattern_prefix(regexstr_prefix);
    const boost::regex pattern_suffix(regexstr_suffix);
 
    int all_prefix_occur(0), all_suffix_occur(0), all_proteins_count(0);
 
    constexpr size_t PROTEIN_CACHE_SIZE = 4e5;
 
    while (true)
    {
      proteins.cacheChunk(PROTEIN_CACHE_SIZE);
      if (!proteins.activateCache()) break;
 
      auto prot_count = (SignedSize)proteins.chunkSize();
      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];
          all_prefix_occur++;
 
          // increase count of observed prefix
          decoy_count[match].first++;
 
          // store observed (case sensitive and with special characters)
          std::string seq_decoy = StringUtils::prefix(seq, match.length());
          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];
          all_suffix_occur++;
 
          // increase count of observed suffix
          decoy_count[match].second++;
 
          // store observed (case sensitive and with special characters)
          std::string seq_decoy = StringUtils::suffix(seq, match.length());
          decoy_case_sensitive[match] = seq_decoy;
        }
      }
    }
 
    // DEBUG ONLY: print counts of found decoys
    for (auto &a : 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 (all_prefix_occur + all_suffix_occur < 0.4 * all_proteins_count)
    {
      OPENMS_LOG_ERROR << "Unable to determine decoy string (not enough occurrences; <40%)!" << std::endl;
      return {false, "?", true};
    }
 
    if (all_prefix_occur == 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_count)
    {
      const std::string & case_insensitive_decoy_string = pair.first;
      const std::pair<int, int>& prefix_suffix_counts = pair.second;
      double freq_prefix = static_cast<double>(prefix_suffix_counts.first) / static_cast<double>(all_prefix_occur);
      double freq_prefix_in_proteins = static_cast<double>(prefix_suffix_counts.first) / static_cast<double>(all_proteins_count);
 
      if (freq_prefix >= 0.8 && freq_prefix_in_proteins >= 0.4)
      {
        if (prefix_suffix_counts.first != 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_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_count)
    {
      const std::string& case_insensitive_decoy_string = pair.first;
      const std::pair<int, int>& prefix_suffix_counts = pair.second;
      double freq_suffix = static_cast<double>(prefix_suffix_counts.second) / static_cast<double>(all_suffix_occur);
      double freq_suffix_in_proteins = static_cast<double>(prefix_suffix_counts.second) / static_cast<double>(all_proteins_count);
 
      if (freq_suffix >= 0.8 && freq_suffix_in_proteins >= 0.4)
      {
        if (prefix_suffix_counts.second != 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_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};
  }
  
private:
  using DecoyStringToAffixCount = std::map<std::string, std::pair<int, int>>;
  using CaseInsensitiveToCaseSensitiveDecoy = std::map<std::string, std::string>;
};
 
} // namespace OpenMS