SequenceContainerTools.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_SEQ_CONTAINER_SEQUENCECONTAINERTOOLS_H
#define BPP_SEQ_CONTAINER_SEQUENCECONTAINERTOOLS_H
 
 
// From the STL:
#include <string>
#include <vector>
#include <map>
#include <memory>
 
#include <Bpp/Numeric/VectorTools.h>
#include "../SymbolListTools.h"
#include "SequenceContainer.h"
#include "VectorSequenceContainer.h"
#include "../Alphabet/CodonAlphabet.h"
 
namespace bpp
{
using SequenceSelection = std::vector<size_t>;
using SiteSelection = std::vector<size_t>;
 
class SequenceContainerTools
{
public:
  SequenceContainerTools() {}
  virtual ~SequenceContainerTools() {}
 
public:
  template<class SequenceType, class HashType>
  static bool hasSequenceWithName(const TemplateSequenceContainerInterface<SequenceType, HashType>& sc, const std::string& name)
  {
    size_t nbSeq = sc.getNumberOfSequences();
    for (size_t i = 0; i < nbSeq; ++i)
    {
      if (sc.sequence(i).getName() == name)
      {
        return true;
      }
    }
    return false;
  }
 
  template<class SequenceType, class HashType>
  static std::unique_ptr< TemplateSequenceContainerInterface<SequenceType, HashType>> createContainerOfSpecifiedSize(std::shared_ptr<const Alphabet>& alphabet, size_t size)
  {
    auto vsc = std::make_unique< TemplateVectorSequenceContainer<SequenceType>>(alphabet);
    for (size_t i = 0; i < size; ++i)
    {
      vsc->addSequence(SequenceType(TextTools::toString(i), "", alphabet), false);
    }
    return vsc;
  }
 
 
  template<class SequenceType, class HashType>
  static std::unique_ptr< TemplateSequenceContainerInterface<SequenceType, HashType>> createContainerWithSequenceNames(
      std::shared_ptr<const Alphabet>& alphabet,
      const std::vector<std::string>& seqNames)
  {
    auto sc = createContainerOfSpecifiedSize<SequenceType, HashType>(alphabet, seqNames.size());
    sc->setSequenceNames(seqNames, true);
    return sc;
  }
 
 
  template<class ContFrom, class ContTo, class SequenceType>
  static void convertContainer(const ContFrom& input, ContTo& output)
  {
    for (size_t i = 0; i < input.getNumberOfSequences(); ++i)
    {
      auto seq = std::make_unique<SequenceType>(input.sequence(i));
      output.addSequence(seq->getName(), seq);
    }
  }
 
  template<class SequenceType, class HashType>
  static void getSelectedSequences(
      const TemplateSequenceContainerInterface<SequenceType, HashType>& sequences,
      const SequenceSelection& selection,
      TemplateSequenceContainerInterface<SequenceType, HashType>& outputCont)
  {
    for (size_t position : selection)
    {
      auto seq = std::make_unique<SequenceType>(sequences.sequence(position));
      outputCont.addSequence(seq->getName(), seq);
    }
  }
 
 
  template<class SequenceType, class HashType>
  static void getSelectedSequences(
      const TemplateSequenceContainerInterface<SequenceType, HashType>& sequences,
      const std::vector<std::string>& selection,
      TemplateSequenceContainerInterface<SequenceType, HashType>& outputCont,
      bool strict = true)
  {
    for (const std::string& key : selection)
    {
      if (strict)
      {
        auto seq = std::make_unique<SequenceType>(sequences.sequence(key));
        outputCont.addSequence(seq->getName(), seq);
      }
      else
      {
        if (sequences.hasSequence(key))
        {
          auto seq = std::make_unique<SequenceType>(sequences.sequence(key));
          outputCont.addSequence(seq->getName(), seq);
        }
      }
    }
  }
 
 
  template<class SequenceType, class HashType>
  static void keepOnlySelectedSequences(
      TemplateSequenceContainerInterface<SequenceType, HashType>& sequences,
      const SequenceSelection& selection)
  {
    std::vector<std::string> keys = sequences.getSequenceKeys();
    std::vector<std::string> selectedKeys = VectorTools::extract<std::string>(keys, selection);
    std::vector<std::string> keysToRemove;
    VectorTools::diff(keys, selectedKeys, keysToRemove);
    for (const std::string& key : keysToRemove)
    {
      // We need to do this because after removal the indices will not be the same!
      // another solution would be to sort decreasingly the indices...
      sequences.removeSequence(key);
    }
  }
 
 
  template<class SequenceType, class HashType>
  static bool sequencesHaveTheSameLength(
      const TemplateSequenceContainerInterface<SequenceType,
      HashType>& sc)
  {
    size_t ns = sc.getNumberOfSequences();
    if (ns <= 1)
      return true;
    size_t length = sc.sequence(0).size();
    for (size_t i = 1; i < ns; ++i)
    {
      if (sc.sequence(i).size() != length)
        return false;
    }
    return true;
  }
 
 
  static void getCounts(const SequenceContainerInterface& sc, std::map<int, unsigned int>& f)
  {
    for (size_t i = 0; i < sc.getNumberOfSequences(); ++i)
    {
      const Sequence& seq = sc.sequence(i);
      for (size_t j = 0; j < seq.size(); ++j)
      {
        f[seq[j]]++;
      }
    }
  }
 
 
  static void getFrequencies(
      const SequenceContainerInterface& sc,
      std::map<int, double>& f,
      double pseudoCount = 0)
  {
    double n = 0;
    for (size_t i = 0; i < sc.getNumberOfSequences(); ++i)
    {
      const Sequence& seq = sc.sequence(i);
      SymbolListTools::getCounts(seq, f, true);
      n += static_cast<double>(seq.size());
    }
 
    if (pseudoCount != 0)
    {
      std::shared_ptr<const Alphabet> pA = sc.getAlphabet();
      for (int i = 0; i < static_cast<int>(pA->getSize()); ++i)
      {
        f[i] += pseudoCount;
      }
      n += pseudoCount * static_cast<double>(pA->getSize());
    }
 
    for (auto& i : f)
    {
      i.second = i.second / n;
    }
  }
 
 
  static void getFrequencies(
      const ProbabilisticSequenceContainerInterface& sc,
      std::map<int, double>& f,
      double pseudoCount = 0)
  {
    double n = 0;
    for (size_t i = 0; i < sc.getNumberOfSequences(); ++i)
    {
      const ProbabilisticSequence& seq = sc.sequence(i);
      SymbolListTools::getCounts(seq, f, true);
      n += static_cast<double>(seq.size());
    }
 
    if (pseudoCount != 0)
    {
      std::shared_ptr<const Alphabet> pA = sc.getAlphabet();
      for (int i = 0; i < static_cast<int>(pA->getSize()); ++i)
      {
        f[i] += pseudoCount;
      }
      n += pseudoCount * static_cast<double>(pA->getSize());
    }
 
    for (auto& i : f)
    {
      i.second = i.second / n;
    }
  }
 
 
  static void getFrequencies(
      const SequenceDataInterface& sc,
      std::map<int, double>& f,
      double pseudoCount = 0)
  {
    try
    {
      getFrequencies(dynamic_cast<const SequenceContainerInterface&>(sc), f, pseudoCount);
      return;
    }
    catch (std::bad_cast&) {}
    try
    {
      getFrequencies(dynamic_cast<const ProbabilisticSequenceContainerInterface&>(sc), f, pseudoCount);
    }
    catch (std::bad_cast&)
    {
      throw Exception("SequenceContainerTools::getFrequencies : unsupported SequenceDataInterface implementation.");
    }
  }
 
  template<class SequenceType, class HashType>
  static void append(
      TemplateSequenceContainerInterface<SequenceType, HashType>& seqCont1,
      const TemplateSequenceContainerInterface<SequenceType, HashType>& seqCont2)
  {
    const auto& keys = seqCont2.getSequenceKeys();
    for (size_t i = 0; i < seqCont2.getNumberOfSequences(); ++i)
    {
      auto tm = std::unique_ptr<SequenceType>(seqCont2.sequence(i).clone());
      seqCont1.addSequence(keys[i], tm);
    }
  }
 
 
  template<class SequenceType, class HashType>
  static void mergeByKey(
      const TemplateSequenceContainerInterface<SequenceType, HashType>& seqCont1,
      const TemplateSequenceContainerInterface<SequenceType, HashType>& seqCont2,
      TemplateSequenceContainerInterface<SequenceType, HashType>& outputCont)
  {
    if (seqCont1.getAlphabet()->getAlphabetType() != seqCont2.getAlphabet()->getAlphabetType())
      throw AlphabetMismatchException("SequenceContainerTools::mergeByKey.", seqCont1.getAlphabet(), seqCont2.getAlphabet());
 
    for (const auto& key: seqCont1.getSequenceKeys())
    {
      auto tmp = std::unique_ptr<SequenceType>(seqCont1.sequence(key).clone());
      tmp->append(seqCont2.sequence(key));
      outputCont.addSequence(key, tmp);
    }
  }
 
  template<class SequenceType, class HashType>
  static void
  convertAlphabet(
      const TemplateSequenceContainerInterface<SequenceType, HashType>& seqCont,
      TemplateSequenceContainerInterface<SequenceType, HashType>& outputCont)
  {
    std::vector<std::string> sequenceKeys = seqCont.getSequenceKeys();
    for (size_t i = 0; i < seqCont.getNumberOfSequences(); ++i)
    {
      std::string seqName = seqCont.sequence(i).getName();
      std::string seqKey = sequenceKeys[i];
      auto alpha = outputCont.getAlphabet();
      auto seq = std::unique_ptr<SequenceType>(new SequenceType(seqName, seqCont.sequence(i).toString(), alpha));
      outputCont.addSequence(seqKey, seq);
    }
  }
 
  template<class SequenceType>
  static std::unique_ptr< TemplateSequenceContainerInterface<SequenceType>>
  getCodonPosition(
      const TemplateSequenceContainerInterface<SequenceType, std::string>& sequences,
      size_t pos)
  {
    auto calpha = std::dynamic_pointer_cast<const CodonAlphabet>(sequences.getAlphabet());
    if (!calpha)
      throw AlphabetException("SequenceContainerTools::getCodonPosition. Input sequences should be of type codon.", sequences.getAlphabet());
    auto newcont = std::make_unique< TemplateVectorSequenceContainer<SequenceType>>(calpha->getNucleicAlphabet());
    for (size_t i = 0; i < sequences.getNumberOfSequences(); ++i)
    {
      const SequenceType& seq = sequences.sequence(i);
      std::vector<int> newseq(seq.size());
      for (size_t j = 0; j < seq.size(); ++j)
      {
        newseq[i] = calpha->getNPosition(seq[i], pos);
      }
      std::shared_ptr<const bpp::Alphabet> na = calpha->getNucleicAlphabet();
      auto s = std::make_unique<SequenceType>(seq.getName(), newseq, seq.getComments(), na);
      newcont->addSequence(sequences.getSequenceKeys()[i], s);
    }
    return newcont;
  }
};
} // end of namespace bpp.
#endif // BPP_SEQ_CONTAINER_SEQUENCECONTAINERTOOLS_H