SequenceApplicationTools.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_SEQ_APP_SEQUENCEAPPLICATIONTOOLS_H
#define BPP_SEQ_APP_SEQUENCEAPPLICATIONTOOLS_H
 
#include <map>
#include <string>
#include <memory>
 
#include "../Alphabet/Alphabet.h"
#include "../AlphabetIndex/AlphabetIndex1.h"
#include "../AlphabetIndex/AlphabetIndex2.h"
#include "../GeneticCode/GeneticCode.h"
#include "../Container/SequenceContainer.h"
#include "../Container/VectorSiteContainer.h"
#include "../Container/SiteContainerTools.h"
#include "../SiteTools.h"
 
namespace bpp
{
class SequenceApplicationTools
{
public:
  SequenceApplicationTools() {}
  virtual ~SequenceApplicationTools() {}
 
public:
  static std::unique_ptr<Alphabet> getAlphabet(
      const std::map<std::string, std::string>& params,
      const std::string& suffix = "",
      bool suffixIsOptional = true,
      bool verbose = true,
      bool allowGeneric = false,
      int warn = 1);
 
  static std::unique_ptr<GeneticCode> getGeneticCode(
      std::shared_ptr<const NucleicAlphabet> alphabet,
      const std::string& description);
 
  static std::unique_ptr<AlphabetIndex1> getAlphabetIndex1(
      std::shared_ptr<const Alphabet> alphabet,
      const std::string& description,
      const std::string& message = "Alphabet measure:",
      bool verbose = true);
 
  static std::unique_ptr<AlphabetIndex1> getAlphabetIndex1(
      std::shared_ptr<const CodonAlphabet> alphabet,
      std::shared_ptr<const GeneticCode> gencode,
      const std::string& description,
      const std::string& message = "Alphabet measure:",
      bool verbose = true);
 
 
  static std::unique_ptr<AlphabetIndex2> getAlphabetIndex2(
      std::shared_ptr<const Alphabet> alphabet,
      const std::string& description,
      const std::string& message = "Alphabet distance:",
      bool verbose = true);
 
  static std::unique_ptr<AlphabetIndex2> getAlphabetIndex2(
      std::shared_ptr<const CodonAlphabet> alphabet,
      std::shared_ptr<const GeneticCode> gencode,
      const std::string& description,
      const std::string& message = "Alphabet distance:",
      bool verbose = true);
 
  static std::unique_ptr<SequenceContainerInterface> getSequenceContainer(
      std::shared_ptr<const Alphabet> alpha,
      const std::map<std::string, std::string>& params,
      const std::string& suffix = "",
      bool suffixIsOptional = true,
      bool verbose = true,
      int warn = 1);
 
  static std::unique_ptr<VectorSiteContainer> getSiteContainer(
      std::shared_ptr<const Alphabet> alpha,
      const std::map<std::string, std::string>& params,
      const std::string& suffix = "",
      bool suffixIsOptional = true,
      bool verbose = true,
      int warn = 1);
 
  static std::unique_ptr<ProbabilisticVectorSiteContainer> getProbabilisticSiteContainer(
      std::shared_ptr<const Alphabet> alpha,
      const std::map<std::string, std::string>& params,
      const std::string& suffix = "",
      bool suffixIsOptional = true,
      bool verbose = true,
      int warn = 1);
 
  static std::map<size_t, std::unique_ptr<VectorSiteContainer>>
  getSiteContainers(
      std::shared_ptr<const Alphabet> alpha,
      const std::map<std::string, std::string>& params,
      const std::string& prefix = "input.",
      const std::string& suffix = "",
      bool suffixIsOptional = true,
      bool verbose = true,
      int warn = 1);
 
  static std::map<size_t, std::unique_ptr<ProbabilisticVectorSiteContainer>>
  getProbabilisticSiteContainers(
      std::shared_ptr<const Alphabet> alpha,
      const std::map<std::string, std::string>& params,
      const std::string& prefix = "input.",
      const std::string& suffix = "",
      bool suffixIsOptional = true,
      bool verbose = true,
      int warn = 1);
 
 
  static void restrictSelectedSequencesByName(
      SequenceContainerInterface& allSequences,
      const std::map<std::string, std::string>& params,
      std::string suffix = "",
      bool suffixIsOptional = true,
      bool verbose = true,
      int warn = 1);
 
 
  template<class SiteType, class SequenceType>
  static std::unique_ptr< TemplateVectorSiteContainer<SiteType, SequenceType>>
  getSitesToAnalyse(
      const TemplateSiteContainerInterface<SiteType, SequenceType, std::string>& allSites,
      const std::map<std::string, std::string>& params,
      std::string suffix = "",
      bool suffixIsOptional = true,
      bool gapAsUnknown = true,
      bool verbose = true,
      int warn = 1)
  {
    // Fully resolved sites, i.e. without jokers and gaps:
    std::unique_ptr< TemplateVectorSiteContainer<SiteType, SequenceType>> sitesToAnalyse;
 
    size_t numSeq = allSites.getNumberOfSequences();
 
    std::string option = ApplicationTools::getStringParameter("input.sequence.sites_to_use", params, "complete", suffix, suffixIsOptional, warn);
    if (verbose)
      ApplicationTools::displayResult("Sites to use", option);
 
    if (option == "all")
    {
      sitesToAnalyse = std::make_unique< TemplateVectorSiteContainer<SiteType, SequenceType>>(allSites);
      size_t nbSites = sitesToAnalyse->getNumberOfSites();
 
      std::string maxGapOption = ApplicationTools::getStringParameter("input.sequence.max_gap_allowed", params, "100%", suffix, suffixIsOptional, warn);
 
      double gapCount = 0;
 
      if (maxGapOption[maxGapOption.size() - 1] == '%')
      {
        double gapFreq = TextTools::toDouble(maxGapOption.substr(0, maxGapOption.size() - 1)) / 100;
        gapCount = gapFreq * (int)numSeq;
      }
      else
        gapCount = TextTools::to<int>(maxGapOption) - NumConstants::TINY();
 
      if (gapCount < static_cast<double>(numSeq) - NumConstants::TINY())
      {
        if (verbose)
          ApplicationTools::displayTask("Remove sites with gaps", true);
        for (size_t i = nbSites; i > 0; i--)
        {
          if (verbose)
            ApplicationTools::displayGauge(nbSites - i, nbSites - 1, '=');
 
          if (static_cast<double>(SiteTools::numberOfGaps(sitesToAnalyse->site(i - 1))) > gapCount)
            sitesToAnalyse->deleteSites(i - 1, 1);
        }
        if (verbose)
          ApplicationTools::displayTaskDone();
      }
 
      std::string maxUnresolvedOption = ApplicationTools::getStringParameter("input.sequence.max_unresolved_allowed", params, "100%", suffix, suffixIsOptional, warn);
 
      double unresCount = 0;
 
      if (maxUnresolvedOption[maxUnresolvedOption.size() - 1] == '%')
      {
        double unresFreq = TextTools::toDouble(maxUnresolvedOption.substr(0, maxUnresolvedOption.size() - 1)) / 100;
        unresCount = unresFreq * (int)numSeq;
      }
      else
        unresCount = TextTools::to<double>(maxUnresolvedOption) - NumConstants::TINY();
 
      nbSites = sitesToAnalyse->getNumberOfSites();
 
      if (unresCount < static_cast<double>(numSeq) - NumConstants::TINY())
      {
        if (verbose)
          ApplicationTools::displayTask("Remove unresolved sites", true);
        for (size_t i = nbSites; i > 0; i--)
        {
          if (verbose)
            ApplicationTools::displayGauge(nbSites - i, nbSites - 1, '=');
 
          if (static_cast<double>(SiteTools::numberOfUnresolved(sitesToAnalyse->site(i - 1))) > unresCount)
            sitesToAnalyse->deleteSites(i - 1, 1);
        }
        if (verbose)
          ApplicationTools::displayTaskDone();
      }
    }
    else if (option == "complete")
    {
      sitesToAnalyse = SiteContainerTools::getCompleteSites(allSites);
      size_t nbSites = sitesToAnalyse->getNumberOfSites();
      if (verbose)
        ApplicationTools::displayResult("Complete sites", TextTools::toString(nbSites));
    }
    else if (option == "nogap")
    {
      sitesToAnalyse = SiteContainerTools::getSitesWithoutGaps(allSites);
      size_t nbSites = sitesToAnalyse->getNumberOfSites();
      if (verbose)
        ApplicationTools::displayResult("Sites without gap", TextTools::toString(nbSites));
    }
    else
    {
      throw Exception("Option '" + option + "' unknown in parameter 'sequence.sites_to_use'.");
    }
 
    auto ca = std::dynamic_pointer_cast<const CodonAlphabet>(sitesToAnalyse->getAlphabet());
    if (ca)
    {
      option = ApplicationTools::getStringParameter("input.sequence.remove_stop_codons", params, "no", suffix, true, warn);
      if ((option != "") && verbose)
        ApplicationTools::displayResult("Remove Stop Codons", option);
 
      if (option == "yes")
      {
        std::string codeDesc = ApplicationTools::getStringParameter("genetic_code", params, "Standard", "", true, warn);
        auto nucAlph = ca->getNucleicAlphabet();
        auto gCode = getGeneticCode(nucAlph, codeDesc);
        SiteContainerTools::removeSitesWithStopCodon(*sitesToAnalyse, *gCode);
      }
    }
 
    if (verbose)
      ApplicationTools::displayResult("Number of sites", sitesToAnalyse->getNumberOfSites());
 
    return sitesToAnalyse;
  }
 
 
  static void writeSequenceFile(
      const SequenceContainerInterface& sequences,
      const std::map<std::string, std::string>& params,
      const std::string& suffix = "",
      bool verbose = true,
      int warn = 1);
 
  static void writeAlignmentFile(
      const SiteContainerInterface& sequences,
      const std::map<std::string, std::string>& params,
      const std::string& suffix = "",
      bool verbose = true,
      int warn = 1);
};
} // end of namespace bpp.
#endif // BPP_SEQ_APP_SEQUENCEAPPLICATIONTOOLS_H