src/sequence.h

#ifndef express_sequence_h
#define express_sequence_h

#include <boost/scoped_array.hpp>
#include <string>
#include "frequencymatrix.h"

inline char ctoi(const char c) {
  switch(c) {
    case 'A':
    case 'a':
       return 0;
    case 'C':
    case 'c':
       return 1;
    case 'G':
    case 'g':
       return 2;
    case 'T':
    case 't':
       return 3;
    default:
       return 0;
  }
}
inline char complement(const char c) {
  return c^3;
}

class Sequence {
 public:
  Sequence() {}
  virtual ~Sequence() {}
  virtual void set(const std::string& seq, bool rev) = 0;
  virtual size_t operator[](const size_t index) const = 0;
  virtual size_t get_ref(const size_t index) const = 0;
  virtual void update_est(const size_t index, const size_t nuc, float mass) = 0;
  virtual void update_obs(const size_t index, const size_t nuc, float mass) = 0;
  virtual void update_exp(const size_t index, const size_t nuc, float mass) = 0;
  virtual float get_prob(const size_t index, const size_t nuc) const = 0;
  virtual float get_obs(const size_t index, const size_t nuc) const = 0;
  virtual float get_exp(const size_t index, const size_t nuc) const = 0;
  virtual bool prob() const = 0;
  virtual size_t length() const = 0;
  virtual bool empty() const = 0;
  virtual void calc_p_vals(std::vector<double>& p_vals) const = 0;
  std::string serialize();
};

class SequenceFwd: public Sequence
{
  boost::scoped_array<const char> _ref_seq;
  FrequencyMatrix<float> _est_seq;
  FrequencyMatrix<float> _obs_seq;
  FrequencyMatrix<float> _exp_seq;
  bool _prob;
  size_t _len;

 public:
  SequenceFwd();
  SequenceFwd(const std::string& seq, bool rev, bool prob=false);
  SequenceFwd(const SequenceFwd& other);
  SequenceFwd& operator=(const SequenceFwd& other);
  // The following methods are documented in the abstract Sequence class.
  void set(const std::string& seq, bool rev);
  size_t operator[](const size_t index) const;
  size_t get_ref(const size_t index) const;
  float get_exp(const size_t index, const size_t nuc) const;
  float get_obs(const size_t index, const size_t nuc) const;
  void update_est(const size_t index, const size_t nuc, float mass);
  void update_obs(const size_t index, const size_t nuc, float mass);
  void update_exp(const size_t index, const size_t nuc, float mass);
  float get_prob(const size_t index, const size_t nuc) const;
  bool prob() const { return _prob; }
  bool empty() const { return _len==0; }
  size_t length() const { return _len; }
  void calc_p_vals(std::vector<double>& p_vals) const;
};

class SequenceRev: public Sequence {
  SequenceFwd* _seq;

 public:
  SequenceRev() : _seq(NULL){}
  SequenceRev(SequenceFwd& seq) : _seq(&seq) {}
  size_t length() const {
    if (_seq == NULL) {
   return 0;
    }
      return _seq->length();
  }
  bool empty() const {
    if (_seq == NULL) {
      return true;
    }
    return _seq->empty();
  }
  // Set is not allowed in this class.
  void set(const std::string& seq, bool rev) { assert(false); exit(1); }
  size_t operator[](const size_t index) const {
    return complement(_seq->operator[](length()-index-1)); }
  size_t get_ref(const size_t index) const {
    return complement(_seq->get_ref(length()-index-1)); }
  float get_obs(const size_t index, const size_t nuc) const {
    return _seq->get_obs(length()-index-1, complement(nuc)); }
  float get_exp(const size_t index, const size_t nuc) const {
    return _seq->get_exp(length()-index-1, complement(nuc)); }
  void update_est(const size_t index, const size_t nuc, float mass) {
    _seq->update_est(length()-index-1, complement(nuc), mass); }
  void update_obs(const size_t index, const size_t nuc, float mass) {
    _seq->update_obs(length()-index-1, complement(nuc), mass); }
  void update_exp(const size_t index, const size_t nuc, float mass) {
    _seq->update_exp(length()-index-1, complement(nuc), mass); }
  float get_prob(const size_t index, const size_t nuc) const {
    return _seq->get_prob(length()-index-1, complement(nuc)); }
  bool prob() const { return _seq->prob(); }
  void calc_p_vals(std::vector<double>& p_vals) const;
};

#endif