SequenceTools.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_SEQ_SEQUENCETOOLS_H
#define BPP_SEQ_SEQUENCETOOLS_H
 
#include <Bpp/Exceptions.h>
#include <Bpp/Numeric/Random/RandomTools.h>
#include <Bpp/Numeric/Stat/StatTest.h>
 
#include "Alphabet/Alphabet.h"
#include "Alphabet/DNA.h"
#include "Alphabet/RNA.h"
#include "Alphabet/RNY.h"
#include "GeneticCode/GeneticCode.h"
#include "NucleicAcidsReplication.h"
#include "Sequence.h"
#include "SymbolListTools.h"
 
// From the STL:
#include <string>
#include <map>
#include <vector>
#include <algorithm>
#include <memory>
 
namespace bpp
{
class BowkerTest :
  public StatTest
{
private:
  double pvalue_;
  double stat_;
 
public:
  BowkerTest() : pvalue_(1.),
    stat_(0.) {}
 
  virtual ~BowkerTest() {}
 
  BowkerTest* clone() const { return new BowkerTest(*this); }
 
public:
  std::string getName() const { return "Bowker's test for homogeneity."; }
  double getStatistic() const { return stat_; }
  double getPValue() const { return pvalue_; }
 
  void setStatistic(double stat) { stat_ = stat; }
  void setPValue(double pvalue) { pvalue_ = pvalue; }
};
 
class SequenceTools :
  public SymbolListTools
{
private:
  static std::shared_ptr<RNY> RNY_;
  static NucleicAcidsReplication DNARep_;
  static NucleicAcidsReplication RNARep_;
  static NucleicAcidsReplication transc_;
 
public:
  SequenceTools() {}
  virtual ~SequenceTools() {}
 
public:
  static bool areSequencesIdentical(const SequenceInterface& seq1, const SequenceInterface& seq2);
 
  static void subseq(const SequenceInterface& sequence, size_t begin, size_t end, SequenceInterface& output)
  {
    if (end < begin || end >= sequence.size())
      throw Exception("SequenceTools::subseq. Invalid coordinates begin=" + TextTools::toString(begin) + ", end=" + TextTools::toString(end) + " for a sequence of size " + TextTools::toString(sequence.size()) + ".");
    std::vector<int> content(end - begin + 1);
    for (size_t i = 0; i <= end - begin; ++i)
    {
      content[i] = sequence[begin + i];
    }
    output.append(content);
  }
 
  template<class SequenceTypeOut>
  static std::unique_ptr<SequenceTypeOut> subseq(const SequenceInterface& sequence, size_t begin, size_t end)
  {
    auto alphaPtr = sequence.getAlphabet();
    auto seq = std::make_unique<SequenceTypeOut>(alphaPtr);
    seq->setName(sequence.getName());
    seq->setComments(sequence.getComments());
    subseq(sequence, begin, end, *seq);
    return seq;
  }
 
 
  template<class SequenceTypeOut>
  static std::unique_ptr<SequenceTypeOut> concatenate(const SequenceInterface& seq1, const SequenceInterface& seq2)
  {
    // Sequence's alphabets matching verification
    if ((seq1.alphabet().getAlphabetType()) != (seq2.alphabet().getAlphabetType()))
      throw AlphabetMismatchException("SequenceTools::concatenate : Sequence's alphabets don't match ", seq1.getAlphabet(), seq2.getAlphabet());
 
    // Sequence's names matching verification
    if (seq1.getName() != seq2.getName())
      throw Exception ("SequenceTools::concatenate : Sequence's names don't match");
 
    // Concatenate sequences and send result
    auto concat = std::make_unique<SequenceTypeOut>(seq1);
    concat->setToSizeR(seq1.size() + seq2.size());
    for (size_t i = 0; i < seq2.size(); ++i)
    {
      (*concat)[seq1.size() + i] = seq2[i];
    }
    return concat;
  }
 
  static void complement(SequenceInterface& seq);
 
  static std::unique_ptr<Sequence> getComplement(const SequenceInterface& sequence);
 
  static std::unique_ptr<Sequence> transcript(const Sequence& sequence);
 
  static std::unique_ptr<Sequence> reverseTranscript(const Sequence& sequence);
 
  static void invert(SequenceInterface& seq);
 
  static std::unique_ptr<SequenceInterface> getInvert(const SequenceInterface& sequence);
 
  static void invertComplement(SequenceInterface& seq);
 
  static double getPercentIdentity(const SequenceInterface& seq1, const SequenceInterface& seq2, bool ignoreGaps = false);
 
  static size_t getNumberOfSites(const SequenceInterface& seq);
 
  static size_t getNumberOfCompleteSites(const SequenceInterface& seq);
 
  static std::unique_ptr<SequenceInterface> getSequenceWithCompleteSites(const SequenceInterface& seq);
 
  static size_t getNumberOfUnresolvedSites(const SequenceInterface& seq);
 
 
  static void removeGaps(SequenceInterface& seq);
 
  static std::unique_ptr<SequenceInterface> getSequenceWithoutGaps(const SequenceInterface& seq);
 
  static void removeStops(SequenceInterface& seq, const GeneticCode& gCode);
 
  static std::unique_ptr<SequenceInterface> getSequenceWithoutStops(const SequenceInterface& seq, const GeneticCode& gCode);
 
  static void replaceStopsWithGaps(SequenceInterface& seq, const GeneticCode& gCode);
 
  static std::unique_ptr<BowkerTest> bowkerTest(const SequenceInterface& seq1, const SequenceInterface& seq2);
 
  static void getPutativeHaplotypes(
      const SequenceInterface& seq,
      std::vector<std::unique_ptr<SequenceInterface>>& hap,
      unsigned int level = 2);
 
  static std::unique_ptr<Sequence> combineSequences(
      const SequenceInterface& s1,
      const SequenceInterface& s2);
 
  static std::unique_ptr<Sequence> subtractHaplotype(
      const SequenceInterface& s,
      const SequenceInterface& h,
      std::string name = "",
      unsigned int level = 1);
 
  static std::unique_ptr<Sequence> RNYslice(const SequenceInterface& sequence, int ph);
 
  static std::unique_ptr<Sequence> RNYslice(const SequenceInterface& sequence);
 
  static void getCDS(SequenceInterface& sequence, const GeneticCode& gCode, bool checkInit, bool checkStop, bool includeInit = true, bool includeStop = true);
 
  static size_t findFirstOf(const SequenceInterface& seq, const SequenceInterface& motif, bool strict = true);
 
  static std::unique_ptr<Sequence> getRandomSequence(std::shared_ptr<const Alphabet>& alphabet, size_t length);
};
} // end of namespace bpp.
#endif // BPP_SEQ_SEQUENCETOOLS_H