bioseq.h

Go to the documentation of this file.

#ifndef BIOSEQ_H
#define BIOSEQ_H

#include "codon.h"
#include <string>
#include <iostream>
#include <map>
#include <algorithm>
#include <list>
#include "Range.h"
#include <vector>
#include "Interval.h"

using namespace std;

class greaterRangeLength {
   public:
      bool operator()(const Range &r1, const Range &r2) const {
         return r1.length() > r2.length();
      }
};

class bioseqexception : public exception {
   public:
      bioseqexception() throw() : errstr("bioseq error") { }
      bioseqexception(const string &msg) throw() : errstr(msg) { }
      ~bioseqexception() throw() { }
      const char* what() const throw() { return errstr.c_str(); }
   private:
      string errstr;
};

// some helper functions
void printFasta(ostream &ous, const string& seq, int width=70);
string reverseComplement(const string& seq);
void reverseComplementInPlace(string& seq);
bool loadFastaIntoMap(const string &file, map<string,string> &store);
enum sequenceType {PROTEINSEQ=1, DNASEQ=2, RNASEQ=4, NUCLEIC_ACID=6, GENERIC=7, UNKNOWN=0};
void translate(string &pep, const string &seq, int begin, int end=-1);
int countInternalStops(const string &seq);
void longestORFPlus(const string &rna, string &pep, int &b, int &e, int &f);
void longestORFPlus(const string &rna, string &pep, int &b, int &e);
void longestORFPlusSuffix(const string &rna, string &pep, int &b, int &e);
void longestORFPlusPrefix(const string &rna, string &pep, int &b, int &e);
bool longestNoStartORFPlus(const string &rna, pair<int,int> &nterm, string &npep,
      pair<int,int> &full, string &pep);
bool longestNoStopORFPlus(const string &rna, pair<int,int> &cterm, string &cpep,
      pair<int,int> &full, string &pep);

//pair<int,int> longestORFIndex(const Protein &p);

Range P2Rindex(const Range &ofb, const int frame);
Interval P2Rindex(const Interval &ofb, const int frame);
int P2Rindex(const int pos, const int frame);
//{ return 3*pos+frame; }

bool overlayRVR(const Range &r, const vector<Range> &vr);
vector<Range> findAllORFIndex(const string &rna, int base=1,
      int pep_lencut=22, int HHcut=160, int HTcut=120, int TTcut=10);
void findAllPepORFIndex(list<Range> &orfrange, const string &ss, int cutoff);


class bioseq {
   public:
      bioseq() : seq(), name(), title(), code(0) {}
      bioseq(const bioseq &s) : seq(s.seq), name(s.name),
                                title(s.title), code(0) { }
      bioseq(const string &s) : seq(s), name(),title(), code(0) {}
      bioseq(const string &n, const string &s) 
         : seq(s), name(n), title(), code(0) {}
      bioseq(const string &n, const string &s, const string &t) 
         : seq(s), name(n), title(t), code(0) {}
      //virtual ~bioseq() { 
      ~bioseq() { if (code != 0) delete[] code; }

      void clear() { if (code != 0) delete[] code; code=0; seq.clear(); name.clear(); title.clear(); }

      void assign(const string &s) { 
         seq=s; 
         if (code != 0) {
            delete[] code;
            code=0;
         } 
      }
      
      void read(const string &file);

      bool read(istream &ins, string &header) throw (bioseqexception);

      string& importFasta(const string &fas);


      const string& toString() const { return seq; }

      void printFastaWithHeader(ostream &ous, int width=70);
      void printFasta(ostream &ous, int width=70);
      friend ostream& operator<<(ostream &ous, const bioseq &s);

      string substr(int b, int e) const; 
      string substring(const int b, const int len) const throw (bioseqexception)
      { if (b>=length() || b < 0) throw bioseqexception("bioseq::substring start index out of range"); }
      string substring(const int b) const throw (bioseqexception)
      { if (b>=length() || b < 0) throw bioseqexception("bioseq::substring start index out of range");
         return seq.substr(b); }

      bioseq subseq(int b, int e) const throw (bioseqexception);

      bioseq subsequence(int b, int len=-1) const throw (bioseqexception);
      void setName(const string &n) { name=n; }
      void setTitle(const string &t) { title=t; }
      const string& getName() const { return name; }
      const string& getTitle() const { return title; }
      const string& getSequence() const { return seq; }

      void rvc();
      map<char, double> getFrequency() const;

      pair<double,double> computeEntropy() const {
         return entropy(seq); }
      pair<double,double> computeEntropy(int b, int e) const throw(bioseqexception) {
         if (b<0 || e>seq.length()-1) { 
            cerr << "index error in computeEntropy()"; 
            throw bioseqexception("index out of range inside bioseq::computeEntropy");
            //exit(1);
         }
         return entropy(seq.substr(b,e-b+1)); }

      static pair<double,double> entropy(const string &seq);
      int length() const { return seq.length(); }
      const int* getcode() const;

      char operator[](const int i) const { return seq[i]; }
      bioseq& operator=(const bioseq &s);
      bioseq& operator+=(const bioseq &s) { seq += s.seq; return *this; }
      bioseq& operator=(const string &str);
      bool empty() { if (seq.empty()) return true; else return false; }
      void randomize() { 
         random_shuffle(seq.begin(), seq.end()); 
         delete[] code;
         code=0; 
      }
      void toLowerCase();
      void toUpperCase();
      virtual sequenceType getSequenceType() const { return GENERIC; }

   protected:
      string seq;
      string name;
      string title;
      mutable int* code;
};

int aachar2num(char a);
char aanum2char(int c);

class Protein : public bioseq {
   public:
      Protein() : bioseq() { codefreq[0]=-1; }
      Protein(const Protein &s) : bioseq(s)  { 
         for (int i=0; i<numcodes; i++) codefreq[i]=s.codefreq[i]; 
      }
      Protein(const string &str) : bioseq(str) { codefreq[0]=-1; }
      Protein(const bioseq &s) : bioseq(s)  { codefreq[0]=-1; }
      Protein(const string &n, const string &s) 
         : bioseq (n,s)  { codefreq[0]=-1; }
      Protein subseq(int b, int e) const {
         return bioseq::subseq(b,e); }
      ~Protein() { /* parent destructor is enough */ }
      Protein& operator=(const Protein &p);
      Protein& operator=(const string &str);

      static const char* threeLetterCode(const char A) 
      { return symbols[2*(toupper(A)-'A')]; }
      static const char* aminoAcidFullName(const char A) 
      { return symbols[2*(toupper(A)-'A')+1]; }

      double relativeEntropy() const;
      double relativeEntropyUniform() const;
      //const int* encode() const;
      //const int* getcode() const { if (code == 0) encode(); return code; }
      const double* getCodeFrequence();
      static const int numcodes=27; // including *
      static const char* oneLetterSymbols;
      sequenceType getSequenceType() const { return PROTEINSEQ; }
      bool hasStart() const { return seq[0] == 'M'; }
      bool hasStop() const { return seq[seq.length()-1] == '*'; }
      bool hasInternalStop() const;
      int countInternalStops();

   private:
      mutable double codefreq[27];
      //mutable int* code;
      static const char* symbols[];
      static const double aafreq[];
      //static const double aafrequniform[];
};

class DNA : public bioseq {
   public:
      DNA() : bioseq() { }
      DNA(const DNA &s) : bioseq(s) { }
      DNA(const string &s) : bioseq(s)  { }
      DNA(const string &n, const string &s) 
         : bioseq(n,s)  {}
      DNA(const string &n, const string &s, const string &t) 
         : bioseq(n,s,t)  { }

      DNA(const bioseq &s) : bioseq(s)  { }
      //~DNA() { if (code != 0) delete[] code; }
      //~DNA() { bioseq::~bioseq(); }
      ~DNA() { /* parent destructor is enought */ }
      DNA& operator=(const DNA &s);
      //DNA& operator+=(const DNA &s) { return bioseq::operator+=(s); }
      DNA& operator=(const string &str);

      Protein translate() const;
      Protein translate(int begin, int end=-1) const;
      DNA subseq(int b, int e) const {
         return bioseq::subseq(b,e); }
      const int* getcodeNuc() const;
      sequenceType getSequenceType() const { return DNASEQ; }
      double GCContent() const;
      double GCContent(long int &A, long int &C, long int &G, long int &T, long int &N) const;
      //void revcomp() { seq=reverseComplement(seq); }
      void revcomp() { reverseComplementInPlace(seq); }
      Protein longestORFForward(int &b, int &e) const;
      void longestORFForward(Protein &p, int &b, int &e) const;
      Protein longestORFReverse(int &b, int &e) const;
      void longestORFReverse(Protein &p, int &b, int &e) const;
      Protein longestORF(char strand, int &b, int &e) const;
      Protein longestORF(int &b, int &e) const;
      void longestORF(Protein &p, int &b, int &e) const;

      static const codon& getCodonTable() { return codontable; }
      static void setCodonTable(const int ctid) { codontable.use(ctid); }

   protected:
      static codon codontable;
      //mutable int* code;
};

class mRNA : public DNA {
   public: 
      mRNA() : DNA() { }
      mRNA(const string &s);

   private:
      int cdsb;
      int cdse;
      Protein prt;
};


#endif