SiteContainerTools.h

Go to the documentation of this file.

#ifndef BPP_SEQ_CONTAINER_SITECONTAINERTOOLS_H
#define BPP_SEQ_CONTAINER_SITECONTAINERTOOLS_H
 
// Copyright or © or Copr. Bio++ Development Team, (November 17, 2004)
// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#include "SiteContainer.h"
#include "VectorSiteContainer.h"
#include "AlignedSequenceContainer.h"
#include "SequenceContainerTools.h"
#include "AlignmentData.h"
#include "../AlphabetIndex/AlphabetIndex2.h"
#include "../DistanceMatrix.h"
#include "../GeneticCode/GeneticCode.h"
#include "../SiteTools.h"
#include "../CodonSiteTools.h"
#include "../Site.h"
#include <Bpp/Numeric/Matrix/Matrix.h>
#include <Bpp/Numeric/Random/RandomTools.h>
 
// From the STL:
#include <vector>
#include <map>
#include <memory>
 
namespace bpp
{
using SiteSelection = std::vector<size_t>;
 
class SiteContainerTools
{
public:
  SiteContainerTools() {}
  virtual ~SiteContainerTools() {}
 
public:
  template<class SiteType, class SequenceType>
  static std::unique_ptr<TemplateVectorSiteContainer<SiteType, SequenceType>>
  getSitesWithoutGaps(const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites)
  {
    std::vector<std::string> sequenceKeys = sites.getSequenceKeys();
    std::shared_ptr<const Alphabet> alphaPtr = sites.getAlphabet();
    auto selectedSites = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(sequenceKeys, alphaPtr);
    for (size_t i = 0; i < sites.getNumberOfSites(); ++i)
    {
      if (!SiteTools::hasGap(sites.site(i)))   // This calls the method dedicated to basic sites
      {
        std::unique_ptr<SiteType> sitePtr(sites.site(i).clone());
        selectedSites->addSite(sitePtr, false);
      }
    }
    selectedSites->setSequenceNames(sites.getSequenceNames(), false);
    return selectedSites;
  }
 
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr<TemplateVectorSiteContainer<SiteType, SequenceType>>
  getCompleteSites(const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites)
  {
    std::vector<std::string> sequenceKeys = sites.getSequenceKeys();
    std::shared_ptr<const Alphabet> alphaPtr = sites.getAlphabet();
    auto selectedSites = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(sequenceKeys, alphaPtr);
    for (size_t i = 0; i < sites.getNumberOfSites(); ++i)
    {
      if (SiteTools::isComplete(sites.site(i)))   // This calls the method dedicated to basic sites
      {
        std::unique_ptr<SiteType> sitePtr(sites.site(i).clone());
        selectedSites->addSite(sitePtr, false);
      }
    }
    selectedSites->setSequenceNames(sites.getSequenceNames(), false);
    return selectedSites;
  }
 
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr<TemplateSiteContainerInterface<SiteType, SequenceType, std::string>>
  removeGapOnlySites(const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites)
  {
    if (sites.getNumberOfSequences() == 0)
      throw Exception("SiteContainerTools::removeGapOnlySites. Container is empty.");
    std::vector<std::string> sequenceKeys = sites.getSequenceKeys();
    auto alphaPtr = sites.getAlphabet();
    auto newContainer = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(sequenceKeys, alphaPtr);
    for (size_t i = 0; i < sites.getNumberOfSites(); ++i)
    {
      const Site& site = sites.site(i);
      if (!SiteTools::isGapOnly(site))
      {
        auto site2 = std::unique_ptr<SiteType>(site.clone());
        newContainer->addSite(site2, false);
      }
    }
    return newContainer;
  }
 
 
  template<class SiteType, class SequenceType, class HashType>
  static void removeGapOnlySites(TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& sites)
  {
    if (sites.getNumberOfSequences() == 0)
      throw Exception("SiteContainerTools::removeGapOnlySites. Container is empty.");
 
    size_t n = sites.getNumberOfSites();
    size_t i = n;
    while (i > 1)
    {
      ApplicationTools::displayGauge(n - i + 1, n);
      const SiteType* site = &sites.site(i - 1); // Note (jdutheil 18/12/22: for some reason a ref here does not work, resorting to pointer)
      if (SiteTools::isGapOnly(*site))
      {
        size_t end = i;
        while (SiteTools::isGapOnly(*site) && i > 1)
        {
          --i;
          site = &sites.site(i - 1);
        }
        sites.deleteSites(i, end - i);
      }
      else
      {
        --i;
      }
    }
    ApplicationTools::displayGauge(n, n);
    const Site& site = sites.site(0);
    if (SiteTools::isGapOnly(site))
      sites.deleteSite(0);
  }
 
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr<TemplateVectorSiteContainer<SiteType, SequenceType>>
  removeGapOrUnresolvedOnlySites(const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites)
  {
    if (sites.getNumberOfSequences() == 0)
      throw Exception("SiteContainerTools::removeGapOrUnresolvedOnlySites. Container is empty.");
 
    std::vector<std::string> sequenceKeys = sites.getSequenceKeys();
    auto alphaPtr = sites.getAlphabet();
    auto newContainer = std::make_unique<TemplateVectorSiteContainer<SiteType, SequenceType>>(sequenceKeys, alphaPtr);
    for (size_t i = 0; i < sites.getNumberOfSites(); ++i)
    {
      const Site& site = sites.site(i);
      if (!SiteTools::isGapOrUnresolvedOnly(site))
      {
        auto site2 = std::unique_ptr<SiteType>(site.clone());
        newContainer->addSite(site2, false);
      }
    }
    return newContainer;
  }
 
 
  template<class SiteType, class SequenceType, class HashType>
  static void removeGapOrUnresolvedOnlySites(TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& sites)
  {
    if (sites.getNumberOfSequences() == 0)
      throw Exception("SiteContainerTools::removeGapOrUnresolvedOnlySites. Container is empty.");
 
    size_t n = sites.getNumberOfSites();
    size_t i = n;
    while (i > 1)
    {
      ApplicationTools::displayGauge(n - i + 1, n);
      const SiteType& site = sites.site(i - 1);
      if (SiteTools::isGapOnly(site))
      {
        size_t end = i;
        while (SiteTools::isGapOrUnresolvedOnly(site) && i > 1)
        {
          --i;
          site = &sites.site(i - 1);
        }
        sites.deleteSites(i, end - i);
      }
      else
      {
        --i;
      }
    }
    ApplicationTools::displayGauge(n, n);
    const SiteType& site = sites.site(0);
    if (SiteTools::isGapOrUnresolvedOnly(site))
      sites.deleteSite(0);
  }
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr< TemplateVectorSiteContainer<SiteType, SequenceType>>
  removeGapSites(
      const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites,
      double maxFreqGaps)
  {
    if (sites.getNumberOfSequences() == 0)
      throw Exception("SiteContainerTools::removeGapSites. Container is empty.");
 
    std::vector<std::string> sequenceKeys = sites.getSequenceKeys();
    auto newContainer = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(sequenceKeys, sites.getAlphabet());
    for (size_t i = 0; i < sites.getNumberOfSites(); ++i)
    {
      std::map<int, double> freq;
      const Site& site = sites.site(i);
      SiteTools::getFrequencies(site, freq);
      if (freq[-1] <= maxFreqGaps)
      {
        auto site2 = std::make_unique<SiteType>(site.clone());
        newContainer->addSite(site2, false);
      }
    }
    newContainer->setSequenceNames(sites.getSequenceNames(), false);
    return newContainer;
  }
 
 
  template<class SiteType, class SequenceType, class HashType>
  static void removeGapSites(
      TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& sites,
      double maxFreqGaps)
  {
    if (sites.getNumberOfSequences() == 0)
      throw Exception("SiteContainerTools::removeGapSites. Container is empty.");
 
    for (size_t i = sites.getNumberOfSites(); i > 0; --i)
    {
      std::map<int, double> freq;
      SiteTools::getFrequencies(sites.site(i - 1), freq);
      if (freq[-1] > maxFreqGaps)
      {
        sites.deleteSite(i - 1);
      }
    }
  }
 
 
  static std::unique_ptr<SiteContainerInterface> getSitesWithoutStopCodon(
      const SiteContainerInterface& sites,
      const GeneticCode& gCode)
  {
    std::shared_ptr<const CodonAlphabet> pca = std::dynamic_pointer_cast<const CodonAlphabet>(sites.getAlphabet());
    if (!pca)
      throw AlphabetException("Not a Codon Alphabet", sites.getAlphabet().get());
    if (sites.getNumberOfSequences() == 0)
      throw Exception("SiteContainerTools::getSitesWithoutStopCodon. Container is empty.");
 
    std::vector<std::string> sequenceKeys = sites.getSequenceKeys();
    auto alphaP = sites.getAlphabet();
    auto newContainer = std::make_unique<VectorSiteContainer>(sequenceKeys, alphaP);
    for (size_t i = 0; i < sites.getNumberOfSites(); ++i)
    {
      const Site& site = sites.site(i);
      if (!CodonSiteTools::hasStop(site, gCode))
      {
        std::unique_ptr<Site> site2(site.clone());
        newContainer->addSite(site2, false);
      }
    }
    newContainer->setSequenceNames(sites.getSequenceNames(), false);
    return newContainer;
  }
 
  static void removeSitesWithStopCodon(
      SiteContainerInterface& sites,
      const GeneticCode& gCode)
  {
    std::shared_ptr<const CodonAlphabet> pca = std::dynamic_pointer_cast<const CodonAlphabet>(sites.getAlphabet());
    if (!pca)
      throw AlphabetException("Not a Codon Alphabet", sites.getAlphabet().get());
    if (sites.getNumberOfSequences() == 0)
      throw Exception("SiteContainerTools::removeSitesWithStopCodon. Container is empty.");
 
    for (size_t i = sites.getNumberOfSites(); i > 0; --i)
    {
      const Site& site = sites.site(i - 1);
      if (CodonSiteTools::hasStop(site, gCode))
        sites.deleteSite(i - 1);
    }
  }
 
  static void removeSitesWithStopCodon(ProbabilisticSiteContainerInterface& sites, const GeneticCode& gCode)
  {
    throw Exception("SiteContainerTools::removeSitesWithStopCodon. Method not supported for probabilistic sequences.");
  }
 
  template<class SiteType, class SequenceType, class HashType>
  static void getSelectedSites(
      const TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& sites,
      const SiteSelection& selection,
      TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& outputSites)
  {
    for (auto pos : selection)
    {
      auto sitePtr = std::unique_ptr<SiteType>(sites.site(pos).clone());
      outputSites.addSite(sitePtr, false);
    }
    outputSites.setSequenceNames(sites.getSequenceNames(), true);
  }
 
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr< TemplateVectorSiteContainer<SiteType, SequenceType>>
  getSelectedSites(
      const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites,
      const SiteSelection& selection)
  {
    auto alphaPtr = sites.getAlphabet();
    auto outputSites = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(sites.getSequenceKeys(), alphaPtr);
    outputSites->setComments(sites.getComments());
    getSelectedSites<SiteType, SequenceType, std::string>(sites, selection, *outputSites);
    return outputSites;
  }
 
 
  static std::unique_ptr<AlignmentDataInterface>
  getSelectedSites(
      const AlignmentDataInterface& sites,
      const SiteSelection& selection)
  {
    try
    {
      auto& sc = dynamic_cast<const SiteContainerInterface&>(sites);
      auto sel = getSelectedSites<Site, Sequence>(sc, selection);
      return std::move(sel);
    }
    catch (std::bad_cast& e) {}
 
    try
    {
      auto& psc = dynamic_cast<const ProbabilisticSiteContainerInterface&>(sites);
      auto sel = getSelectedSites<ProbabilisticSite, ProbabilisticSequence>(psc, selection);
      return std::move(sel);
    }
    catch (std::bad_cast& e) {}
 
    throw Exception("SiteContainerTools::getSelectedSites : unsupported container type.");
  }
 
 
  template<class SiteType, class SequenceType, class HashType>
  static void getSelectedPositions(
      const TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& sites,
      const SiteSelection& selection,
      TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& outputSites)
  {
    size_t wsize = sites.getAlphabet()->getStateCodingSize();
    if (wsize > 1)
    {
      if (selection.size() % wsize != 0)
        throw IOException("SiteContainerTools::getSelectedPositions: Positions selection is not compatible with the alphabet in use in the container.");
      SiteSelection selection2;
      for (size_t i = 0; i < selection.size(); i += wsize)
      {
        if (selection[i] % wsize != 0)
          throw IOException("SiteContainerTools::getSelectedPositions: Positions selection is not compatible with the alphabet in use in the container.");
 
        for (size_t j = 1; j < wsize; ++j)
        {
          if (selection[i + j] != (selection[i + j - 1] + 1))
            throw IOException("SiteContainerTools::getSelectedPositions: Positions selection is not compatible with the alphabet in use in the container.");
        }
        selection2.push_back(selection[i] / wsize);
      }
      getSelectedSites(sites, selection2, outputSites);
    }
    else
    {
      getSelectedSites(sites, selection, outputSites);
    }
  }
 
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr< TemplateVectorSiteContainer<SiteType, SequenceType>>
  getSelectedPositions(
      const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites,
      const SiteSelection& selection)
  {
    auto alphaPtr = sites.getAlphabet();
    auto outputSites = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(sites.getSequenceKeys(), alphaPtr);
    outputSites->setComments(sites.getComments());
    getSelectedPositions<SiteType, SequenceType, std::string>(sites, selection, *outputSites);
    return outputSites;
  }
 
 
  static std::unique_ptr<Sequence> getConsensus(
      const SiteContainerInterface& sc,
      const std::string& name = "consensus",
      bool ignoreGap = true,
      bool resolveUnknown = false);
 
  static void changeGapsToUnknownCharacters(SiteContainerInterface& sites);
 
  static void changeGapsToUnknownCharacters(ProbabilisticSiteContainerInterface& sites);
 
  static void changeUnresolvedCharactersToGaps(SiteContainerInterface& sites);
 
  static std::unique_ptr<SiteContainerInterface> resolveDottedAlignment(
      const SiteContainerInterface& dottedAln,
      std::shared_ptr<const Alphabet>& resolvedAlphabet);
 
  static std::map<size_t, size_t> getSequencePositions(const Sequence& seq);
 
  static std::map<size_t, size_t> getAlignmentPositions(const Sequence& seq);
 
  static void getSequencePositions(
      const SiteContainerInterface& sites,
      Matrix<size_t>& positions);
  static std::map<size_t, size_t> translateAlignment(
      const Sequence& seq1,
      const Sequence& seq2);
 
  static std::map<size_t, size_t> translateSequence(
      const SiteContainerInterface& sequences,
      size_t i1,
      size_t i2);
 
  static std::unique_ptr<AlignedSequenceContainer> alignNW(
      const Sequence& seq1,
      const Sequence& seq2,
      const AlphabetIndex2& s,
      double gap);
 
  static std::unique_ptr<AlignedSequenceContainer> alignNW(
      const Sequence& seq1,
      const Sequence& seq2,
      const AlphabetIndex2& s,
      double opening,
      double extending);
 
  template<class SiteType, class SequenceType, class HashType>
  static void sampleSites(
      const TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& sites,
      size_t nbSites,
      TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& outSites,
      std::shared_ptr< std::vector<size_t>> index = nullptr)
  {
    for (size_t i = 0; i < nbSites; ++i)
    {
      size_t pos = static_cast<size_t>(RandomTools::giveIntRandomNumberBetweenZeroAndEntry(static_cast<int>(sites.getNumberOfSites())));
      auto s = std::unique_ptr<SiteType>(sites.site(pos).clone());
      outSites.addSite(s, false);
 
      if (index)
        index->push_back(pos);
    }
    outSites.setSequenceNames(sites.getSequenceNames(),true);
  }
 
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr< TemplateVectorSiteContainer<SiteType, SequenceType>>
  sampleSites(
      const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites,
      size_t nbSites,
      std::shared_ptr< std::vector<size_t>> index = nullptr)
  {
    auto sampledSites = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(sites.getAlphabet());
    sampleSites<SiteType, SequenceType, std::string>(sites, nbSites, *sampledSites, index);
    return sampledSites;
  }
 
 
  template<class SiteType, class SequenceType, class HashType>
  static void bootstrapSites(
      const TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& sites,
      TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& outputSites)
  {
    sampleSites(sites, sites.getNumberOfSites(), outputSites, nullptr);
  }
 
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr< TemplateVectorSiteContainer<SiteType, SequenceType>>
  bootstrapSites(const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& sites)
  {
    auto outputSites = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(sites.getAlphabet());
    bootstrapSites<SiteType, SequenceType, std::string>(sites, *outputSites);
    return outputSites;
  }
 
 
  static double computeSimilarity(
      const SequenceInterface& seq1,
      const SequenceInterface& seq2,
      bool dist = false,
      const std::string& gapOption = SIMILARITY_NODOUBLEGAP,
      bool unresolvedAsGap = true);
 
  static std::unique_ptr<DistanceMatrix> computeSimilarityMatrix(
      const SiteContainerInterface& sites,
      bool dist = false,
      const std::string& gapOption = SIMILARITY_NOFULLGAP,
      bool unresolvedAsGap = true);
 
  static const std::string SIMILARITY_ALL;
  static const std::string SIMILARITY_NOFULLGAP;
  static const std::string SIMILARITY_NODOUBLEGAP;
  static const std::string SIMILARITY_NOGAP;
 
  template<class SiteType, class SequenceType, class HashType>
  static void merge(
      TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& seqCont1,
      const TemplateSiteContainerInterface<SiteType, SequenceType, HashType>& seqCont2,
      bool leavePositionAsIs = false)
  {
    if (seqCont1.getAlphabet()->getAlphabetType() != seqCont2.getAlphabet()->getAlphabetType())
      throw AlphabetMismatchException("SiteContainerTools::merge.", seqCont1.getAlphabet(), seqCont2.getAlphabet());
 
    std::vector<HashType> seqKeys1 = seqCont1.getSequenceKeys();
    std::vector<HashType> seqKeys2 = seqCont2.getSequenceKeys();
    const TemplateSiteContainerInterface<SiteType, SequenceType, HashType>* seqCont2bis = 0;
    bool del = false;
    if (seqKeys1 == seqKeys2)
    {
      seqCont2bis = &seqCont2;
    }
    else
    {
      // We shall reorder sequences first:
      auto seqCont2ter = seqCont2.createEmptyContainer();
      SequenceContainerTools::getSelectedSequences(seqCont2, seqKeys1, *seqCont2ter);
      seqCont2bis = seqCont2ter;
      del = true;
    }
 
    if (leavePositionAsIs)
    {
      for (size_t i = 0; i < seqCont2bis->getNumberOfSites(); ++i)
      {
        std::unique_ptr<Site> site(seqCont2bis->site(i).clone());
        seqCont1.addSite(site, false);
      }
    }
    else
    {
      int offset = static_cast<int>(seqCont1.getNumberOfSites());
      for (size_t i = 0; i < seqCont2bis->getNumberOfSites(); ++i)
      {
        std::unique_ptr<Site> site(seqCont2bis->site(i).clone());
        site->setCoordinate(offset + site->getCoordinate());
        seqCont1.addSite(site, false);
      }
    }
 
    if (del)
      delete seqCont2bis;
  }
 
 
  static std::vector<int> getColumnScores(const Matrix<size_t>& positions1, const Matrix<size_t>& positions2, int na = 0);
 
  static std::vector<double> getSumOfPairsScores(const Matrix<size_t>& positions1, const Matrix<size_t>& positions2, double na = 0);
};
} // end of namespace bpp.
#endif // BPP_SEQ_CONTAINER_SITECONTAINERTOOLS_H