CategorySubstitutionRegister.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_MAPPING_CATEGORYSUBSTITUTIONREGISTER_H
#define BPP_PHYL_MAPPING_CATEGORYSUBSTITUTIONREGISTER_H
 
 
#include "SubstitutionRegister.h"
 
namespace bpp
{
class CategorySubstitutionRegister :
  public AbstractSubstitutionRegister
{
protected:
  bool within_;
  size_t nbCategories_;
  mutable std::map<size_t, size_t> categories_;
  std::vector<std::string> categoryNames_;
  std::vector< std::vector<size_t>> index_;
  std::vector< std::vector<size_t>> revIndex_;
 
  bool stationarity_;
 
public:
  CategorySubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap, bool within = false) :
    AbstractSubstitutionRegister(stateMap, "Category"),
    within_(within),
    nbCategories_(0),
    categories_(),
    categoryNames_(),
    index_(),
    revIndex_(),
    stationarity_()
  {}
 
protected:
  template<class T>
  void setAlphabetCategories(const std::map<int, T>& categories)
  {
    // We need to convert alphabet states into model states.
    std::map<size_t, T> modelCategories;
    for (typename std::map<int, T>::const_iterator it = categories.begin(); it != categories.end(); ++it)
    {
      std::vector<size_t> states = stateMap().getModelStates(it->first);
      for (auto st : states)
      {
        modelCategories[st] = it->second;
      }
    }
    // Then we forward the model categories:
    setModelCategories<T>(modelCategories);
  }
 
  template<class T>
  void setModelCategories(const std::map<size_t, T>& categories)
  {
    // First index categories:
    nbCategories_ = 0;
    std::map<T, size_t> cats;
    for (auto it = categories.begin(); it != categories.end(); ++it)
    {
      if (cats.find(it->second) == cats.end())
      {
        ++nbCategories_;
        cats[it->second] = nbCategories_;
      }
    }
 
    // Now creates categories:
    categories_.clear();
    categoryNames_.resize(nbCategories_);
    for (size_t i = 0; i < stateMap().getNumberOfModelStates(); ++i)
    {
      auto it = categories.find(i);
      if (it != categories.end())
      {
        categories_[i] = cats[it->second];
        categoryNames_[cats[it->second] - 1] += stateMap().getStateDescription(i);
      }
      else
      {
        categories_[i] = 0;
      }
    }
 
    size_t count = 1;
    index_.resize(nbCategories_);
    for (size_t i = 0; i < index_.size(); ++i)
    {
      index_[i].resize(nbCategories_);
      for (size_t j = 0; j < index_.size(); ++j)
      {
        if (j != i)
        {
          index_[i][j] = count++;
          std::vector<size_t> pos(2);
          pos[0] = i;
          pos[1] = j;
          revIndex_.push_back(pos);
        }
      }
    }
    if (within_)
    {
      for (size_t i = 0; i < index_.size(); ++i)
      {
        index_[i][i] = count++;
        std::vector<size_t> pos(2);
        pos[0] = i;
        pos[1] = i;
        revIndex_.push_back(pos);
      }
    }
  }
 
public:
  virtual size_t getCategory(size_t state) const
  {
    return categories_[state];
  }
 
  virtual size_t getCategoryFrom(size_t type) const
  {
    if (type <= nbCategories_ * (nbCategories_ - 1))
    {
      return revIndex_[type - 1][0] + 1;
    }
    else
    {
      if (within_)
        return revIndex_[type - 1][0] + 1;
      else
        throw Exception("CategorySubstitutionRegister::getCategoryFrom. Bad substitution type.");
    }
  }
 
  virtual size_t getCategoryTo(size_t type) const
  {
    if (type <= nbCategories_ * (nbCategories_ - 1))
    {
      return revIndex_[type - 1][1] + 1;
    }
    else
    {
      if (within_)
        return revIndex_[type - 1][1] + 1;
      else
        throw Exception("CategorySubstitutionRegister::getCategoryTo. Bad substitution type.");
    }
  }
 
  virtual std::string getCategoryName(size_t category) const
  {
    return categoryNames_[category - 1];
  }
 
  virtual bool allowWithin() const { return within_; }
 
  bool isStationary() const
  {
    return stationarity_;
  }
 
  void setStationarity(bool stat)
  {
    stationarity_ = stat;
  }
 
  size_t getNumberOfCategories() const { return nbCategories_; }
 
  size_t getNumberOfSubstitutionTypes() const { return nbCategories_ * (nbCategories_ - 1) + (within_ ? nbCategories_ : 0); }
 
  virtual size_t getType(size_t fromState, size_t toState) const
  {
    size_t fromCat = categories_[fromState];
    size_t toCat   = categories_[toState];
    if (fromCat > 0 && toCat > 0)
      return index_[fromCat - 1][toCat - 1];
    else
      return 0;
  }
 
  std::string getTypeName(size_t type) const
  {
    return getCategoryName(getCategoryFrom(type)) + "->" +  getCategoryName(getCategoryTo(type));
  }
};
 
 
class ComprehensiveSubstitutionRegister :
  public CategorySubstitutionRegister
{
public:
  ComprehensiveSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap, bool within = false) :
    CategorySubstitutionRegister(stateMap, within)
  {
    std::map<int, int> categories;
    for (int i = 0; i < static_cast<int>(stateMap->getNumberOfModelStates()); ++i)
    {
      categories[i] = i;
    }
    setAlphabetCategories<int>(categories);
  }
 
  ComprehensiveSubstitutionRegister* clone() const { return new ComprehensiveSubstitutionRegister(*this); }
};
 
class GCSubstitutionRegister :
  public CategorySubstitutionRegister
{
public:
  GCSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap, bool within = false) :
    CategorySubstitutionRegister(stateMap, within)
  {
    std::map<int, int> categories;
    categories[0] = 1;
    categories[1] = 2;
    categories[2] = 2;
    categories[3] = 1;
    setAlphabetCategories<int>(categories);
  }
 
  GCSubstitutionRegister* clone() const { return new GCSubstitutionRegister(*this); }
};
 
 
class GCSynonymousSubstitutionRegister :
  public CategorySubstitutionRegister
{
private:
  std::shared_ptr<const GeneticCode> genCode_;
 
public:
  GCSynonymousSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap, std::shared_ptr<const GeneticCode> genCode, bool within = false) :
    CategorySubstitutionRegister(stateMap, within),
    genCode_(genCode)
  {
    auto pCA = genCode_->getCodonAlphabet();
 
    std::map<int, int> categories;
    for (int i = 0; i < static_cast<int>(pCA->getSize()); i++)
    {
      int n = pCA->getThirdPosition(i);
      switch (n)
      {
      case 0:
      case 3:
        categories[i] = 1;
        break;
      case 1:
      case 2:
        categories[i] = 2;
        break;
      }
    }
    setAlphabetCategories<int>(categories);
  }
 
  GCSynonymousSubstitutionRegister(const GCSynonymousSubstitutionRegister& reg) :
    CategorySubstitutionRegister(reg),
    genCode_(reg.genCode_)
  {}
 
  GCSynonymousSubstitutionRegister& operator=(const GCSynonymousSubstitutionRegister& reg)
  {
    CategorySubstitutionRegister::operator=(reg);
    genCode_ = reg.genCode_;
    return *this;
  }
 
  GCSynonymousSubstitutionRegister* clone() const { return new GCSynonymousSubstitutionRegister(*this); }
 
public:
  size_t getType(size_t fromState, size_t toState) const
  {
    int x = stateMap().getAlphabetStateAsInt(fromState);
    int y = stateMap().getAlphabetStateAsInt(toState);
 
    auto pCA = genCode_->getCodonAlphabet();
 
    if (genCode_->isStop(x) || genCode_->isStop( y) || !genCode_->areSynonymous(x, y))
      return 0;
 
    // only substitutions between 3rd positions
 
    if ((pCA->getFirstPosition(x) != pCA->getFirstPosition(y)) ||
        (pCA->getSecondPosition(x) != pCA->getSecondPosition(y)))
      return 0;
 
    size_t fromCat = categories_[fromState];
    size_t toCat   = categories_[toState];
 
    if (fromCat > 0 && toCat > 0)
      return index_[fromCat - 1][toCat - 1];
    else
      return 0;
  }
 
  std::string getTypeName (size_t type) const
  {
    switch (type)
    {
    case 0:
      return "no AT<->GC substitution or non-synonymous substitution";
    case 1:
      return "AT->GC synonymous";
    case 2:
      return "GC->AT synonymous";
    case 3:
      return "AT->AT synonymous";
    case 4:
      return "GC->GC synonymous";
    default:
      throw Exception("GCSynonymousSubstitutionRegister::getTypeName. Bad substitution type.");
    }
  }
};
 
 
class GCPositionSubstitutionRegister :
  public CategorySubstitutionRegister
{
private:
  std::shared_ptr<const GeneticCode> genCode_;
  size_t pos_;
 
public:
  GCPositionSubstitutionRegister(
      std::shared_ptr<const StateMapInterface> stateMap,
      std::shared_ptr<const GeneticCode> genCode,
      size_t pos,
      bool within = false) :
    CategorySubstitutionRegister(stateMap, within),
    genCode_(genCode),
    pos_(pos)
  {
    auto pCA = genCode_->getCodonAlphabet();
 
    std::map<int, int> categories;
    for (int i = 0; i < static_cast<int>(pCA->getSize()); i++)
    {
      int n = pCA->getNPosition(i, pos_);
      switch (n)
      {
      case 0:
      case 3:
        categories[i] = 1;
        break;
      case 1:
      case 2:
        categories[i] = 2;
        break;
      }
    }
    setAlphabetCategories<int>(categories);
  }
 
  GCPositionSubstitutionRegister(const GCPositionSubstitutionRegister& reg) :
    CategorySubstitutionRegister(reg),
    genCode_(reg.genCode_),
    pos_(reg.pos_)
  {}
 
  GCPositionSubstitutionRegister& operator=(const GCPositionSubstitutionRegister& reg)
  {
    CategorySubstitutionRegister::operator=(reg);
    genCode_ = reg.genCode_;
    pos_ = reg.pos_;
    return *this;
  }
 
  GCPositionSubstitutionRegister* clone() const { return new GCPositionSubstitutionRegister(*this); }
 
public:
  size_t getType(size_t fromState, size_t toState) const
  {
    int x = stateMap().getAlphabetStateAsInt(fromState);
    int y = stateMap().getAlphabetStateAsInt(toState);
 
    auto pCA = genCode_->getCodonAlphabet();
 
    if (genCode_->isStop(x) || genCode_->isStop(y))
      return 0;
 
    // only substitutions between Nth positions
 
    for (size_t pos = 0; pos < 3; pos++)
    {
      if (pos != pos_)
        if (pCA->getNPosition(x, pos) != pCA->getNPosition(y, pos))
          return 0;
    }
 
    size_t fromCat = categories_[fromState];
    size_t toCat   = categories_[toState];
 
    if (fromCat > 0 && toCat > 0)
      return index_[fromCat - 1][toCat - 1];
    else
      return 0;
  }
 
  std::string getTypeName (size_t type) const
  {
    auto p(TextTools::toString(pos_ + 1));
 
    switch (type)
    {
    case 0:
      return "no AT" + p + "<->GC" + p + " substitution";
    case 1:
      return "AT" + p + "->GC" + p;
    case 2:
      return "GC" + p + "->AT" + p;
    case 3:
      return "AT" + p + "->AT" + p;
    case 4:
      return "GC" + p + "->GC" + p;
    default:
      throw Exception("GCPositionSubstitutionRegister::getTypeName. Bad substitution type.");
    }
  }
};
} // end of namespace bpp.
#endif // BPP_PHYL_MAPPING_CATEGORYSUBSTITUTIONREGISTER_H