GenModel.h

Go to the documentation of this file.

#ifndef GENMODEL_H
#define GENMODEL_H

#include <string>
#include <vector>
#include "Range.h"
#include <iostream>
#include <cassert>
#include "bioseq.h"
#include <map>
#include <set>

using namespace std;

class PointOutChain : public exception {
   private:
      string message;

   public:
      PointOutChain() throw() : message("Point outside of Chain") { }
      PointOutChain(const char* msg) throw() : message(msg) { }
      explicit PointOutChain(const string &msg) throw() : message(msg) { }
      ~PointOutChain() throw() { }
      virtual const char* what() const throw() { 
         return message.c_str(); }
};

class Badinput : public exception {
   private:
      string message;

   public:
      Badinput() throw() : message("Coordinates bad") { }
      Badinput(const char msg[]) throw() : message(msg) { }
      explicit Badinput(const string &msg) throw() : message(msg) { }
      ~Badinput() throw() { }
      virtual const char* what() const throw() {
         return message.c_str();
      }
};

class InvalidModel : public exception {
   private: 
      string message;
   public:
      InvalidModel() throw() : message("Invalid Model") { }
      InvalidModel(const char msg[]) throw() : message(msg) { }
      ~InvalidModel() throw() { }
      virtual const char* what() const throw() {
         return message.c_str();
      }
};

// helper functions ////////
void string2Vector(vector<int> &v, const string &str);
bool nonincreasingVector(const vector<int> &v);
bool nondecreasingVector(const vector<int> &v);
bool allLessVector(const vector<int> &v1, const vector<int> &v2);

class Noschain : public Range {
   public:
      Noschain() : Range(), exons() { }
      Noschain(const string &exstarts, const string &exends, char strand) throw(Badinput);
      Noschain(const string &exstr) throw(Badinput);
      //Noschain(const string &exstarts, const string &exends, unsigned char strand);
      Noschain(const Noschain &eg);
      /* Assume strand is +;  only match string such as 38M 
       * Cannot handle Hard clipping. Never seen this yet. */
      //Noschain(int begin, const string &cigar);
      /* strand is the direction of the alignment in mRNA direction. */
      //Noschain(int begin, const string &cigar, char strand);
      void fixAndBuild(const vector<int> &bb, const vector<int> &ee, char strand);
      ~Noschain() { }
      Noschain& operator=(const Noschain &eg);
      Noschain& operator+=(const Range &r); 
      Noschain operator+(int c) const;
      Noschain operator-(int c) const;
      friend Noschain operator+(int c, const Noschain& chain) {
         return chain+c; }
      friend Noschain operator-(int c, const Noschain& chain);
      Noschain& operator-=(int c);
      Noschain& operator+=(int c);

      Range& operator[](int idx) { return exons[idx]; }
      const Range& operator[](int idx) const { return exons[idx]; }
      Range firstExon() const { return exons[0]; }
      Range lastExon() const { return exons.back(); }

      bool operator==(const Noschain &nc) const { return exons == nc.exons; }
      bool operator!=(const Noschain &nc) const { return exons != nc.exons; }
      bool operator<(const Noschain &nc) const;
      int compareByDirection(const Noschain &nc) const;
      Noschain commonExons(const Noschain &eg) const;
      Noschain commonIntrons(const Noschain &nc) const;

      bool isNull() const { return Range::isNull() && exons.empty(); }
      bool empty() const { return exons.empty(); }
      void clear() { exons.clear(); Range::setNull(); }
      void setNull() { Range::setNull(); exons.clear(); }
      int exonLength() const;
      int intronLength() const;
      Range maxIntron() const;
      int maxIntronLength() const { return maxIntron().length(); }
      Range minIntron() const;
      int minIntronLength() const { return minIntron().length(); }
      vector<Range> introns() const;
      bool containIntron(int p1, int p2) const {
         int i1=exonIndex(p1); int i2=exonIndex(p2);
         return i1 != -1 && i2 != -1 && i1  != i2; 
      }
      Range outerRange() const {
         //return Range(exons.front().begin(), exons.back().end());
         return Range(b,e);
      }
      vector<Range> intronsInside(const int len) const;

      int numberOfRanges() const { return exons.size(); }
      Noschain subchainBeforePoint(int p) const throw (PointOutChain);
      Noschain subchainAfterPoint(int p) const throw (PointOutChain);
      void trimAfterPoint(const int p) throw (PointOutChain);
      void trimBeforePoint(const int p) throw (PointOutChain);
      Noschain subchain(const Range &r) const throw (PointOutChain);
      Noschain subchain(int bg, int ed) const throw (PointOutChain)
      { return subchain(Range(bg,ed)); }
      const vector<Range>& getExons() const { return exons; }
      void setExons(const vector<Range>& ex) { exons=ex; updateOuterRange(); }
      void setChain(const Noschain &chain) { setExons(chain.exons); }
      void updateOuterRange() { Range::set(exons[0].begin(), exons.back().end()); }
      void set(const Range &r) throw (Badinput);
      void setEnd(int ne) throw (PointOutChain);
      void setBegin(int nb);
      void growEnd(int icr);
      int advancePosOnExon(int pos, unsigned int d) const;
      int retreatPosOnExon(int pos, unsigned int d) const;
      int exonIndex(int p) const;
      int exonIndex(const Range &r) const;
      bool insideExon(int p) const { return exonIndex(p) != -1; }
      bool insideExon(const Range &r, bool samedir=false) const;
      
      /*** Output functions *****/

      virtual ostream& show(ostream &ous=cout) const;
      friend ostream& operator<<(ostream &ous, const Noschain &nc) 
         { return nc.print(ous); }
      virtual ostream& print(ostream &ous) const;
      string toString() const;
      pair<string,string> startEnd(char sep=',') const;
      string jgiformat(char sep='\t') const;
      string sfExonsGenomic() const;
      string sfExonsTranscript() const;

      Noschain onewayChain() const;
      void onewayChain(Noschain &owc) const;
      void exonOnlyChain(Noschain &eoc) const;
      Noschain exonOnlyChain() const;
      string asDelimitedString(const char sep=':', const char rsep[]=",") const;
      string asDelimitedString(const char sep=':', const char rsep=',') const;

      //bool exonContain(const Noschain &nc, bool samedirection=false) const throw(PointOutChain);
      bool exonContain(const Noschain &nc, bool samedirection=false) const;
      int compareChainAndFix(Noschain &nc, int &edit);
      bool exonOverlap(const Noschain &nc) const;
      int exonOverlapLength(const Noschain &nc) const;
      int exonIntersectLength(const Noschain &nc) const;
      pair<float, float> exonIntersectFraction(const Noschain &nc) const;
      bool append(Noschain &nc);
      void extendEnd(const Noschain &nc);

      Noschain& reverse();
      void subsequence(string &sub, const string &gs) const throw(Badinput);
      bool isDangle(const int intronlen=900, const int exonlen=27) const;
      bool directionAgree() const;

   protected:
      vector<Range> exons;
      //Range cds;
};

void errorMessageFix1(Noschain &c1, Noschain &c2, int status, int end);

extern int FUZZYMARGIN; //=16;
int comparePlusChainFirstBefore(Noschain &c1, Noschain &c2);
int comparePlusChainIdenticalLeft(Noschain &c1, Noschain &c2);
int compareMinusChainFirstBefore(Noschain &c1, Noschain &c2);
int compareMinusChainMinusIdenticalRight(Noschain &c1, Noschain &c2);

class lessByChainDirectionPtr {
   public:
      bool operator()(const Noschain *c1, const Noschain *c2) const {
         return c1->compareByDirection(*c2) == -1;
      }
};

class lessByChainDirection {
   public:
      bool operator()(const Noschain &c1, const Noschain &c2) const {
         return c1.compareByDirection(c2) == -1;
      }
};

class lessChainPtr {
   public:
      bool operator()(const Noschain *c1, const Noschain *c2) const {
         return *c1 < *c2;
      }
};

class Alnchain : public Noschain {
   public:
      Alnchain(const string &exstr) 
         : Noschain(exstr), iden(100), cov(1) { }
      Alnchain(int begin, const string &cigar, int nmismatch);
      Alnchain(int begin, const string &cigar, char strand, int nmismatch);
      Alnchain() : Noschain(), iden(100), cov(1) { }
      Alnchain(double id, double cv) : Noschain(), iden(id), cov(cv) { }
      Alnchain(const string &exons, double identity, double coverage) : Noschain(exons), iden(identity), cov(coverage) {}

      double identity() const { return iden; }
      double coverage() const { return cov; }
      void setIdentity(double newiden) { iden=newiden; }
      void setCoverage(double newcov) { cov=newcov; }
      string toString(const char sep='\t') const;
      Alnchain& operator=(const Alnchain &ch);
      ostream& print(ostream &ous) const;
      //friend ostream& operator<<(ostream &ous, const Alnchain &chn);

   protected:
      double iden;
      double cov;
};

class Alnchainid : public Alnchain {
   public:
      Alnchainid() : Alnchain(), estids() { }
      Alnchainid(double seqiden, double cv, const string &estid) 
         : Alnchain(seqiden, cv), estids(1, estid) { }
      Alnchainid(int begin, const string &cigar, int nmismatch, const string &estid)
         : Alnchain(begin,cigar,nmismatch), estids(1,estid) { }
      Alnchainid(int begin, const string &cigar, char strand, int nmismatch, const string &estid)
         : Alnchain(begin,cigar,strand,nmismatch), estids(1,estid) { }
      Alnchainid(const string &exons, double identity, double coverage,
              const string &estidstr) : Alnchain(exons, identity, coverage) { setESTIds(estidstr); }

      void addESTId(const string &estid) { estids.push_back(estid); }
      void setESTIds(const string &idstr);
      int numEST() const { return estids.size(); }
      const vector<string>& getESTIds() const { return estids; }
      vector<string>& getESTIds() { return estids; }
      string toString(const char sep='\t') const;
      string ESTIdString(const char sep=',') const;
      Alnchainid& operator=(const Alnchainid &ch);
      void absorbESTIds(const Alnchainid &ch) { estids.insert(estids.end(), ch.estids.begin(), ch.estids.end()); }
      void averageWith(const Alnchainid &chn);
      ostream& print(ostream &ous) const;

   private:
      vector<string> estids;
};


void readandstoreGmap(const string &inputFile,  
      map<string, set<Alnchainid*, lessChainPtr> > &alnseg,
      map<int,int> &pathstat);
void readandstoreGmap(istream &ins,  map<string, set<Alnchainid*, lessChainPtr> > &alnseg,
      map<int,int> &pathstat);

void archiveAln(map<string, set<Alnchainid*, lessChainPtr> > &alnseg, const string &file);
void archiveAln(map<string, set<Alnchainid*, lessChainPtr> > &alnseg, ostream &ous);

class GenModel : public Noschain {
   public:
      GenModel() : Noschain(), cdsphase(0), numstop(0),
               hasstart(true), hasstop(true),
               pep(0) { }
      GenModel(const string &exstarts, const string &exends, 
            char strand, int cdsb, int cdse, const string &gi, 
            int ii, const string &si) throw(InvalidModel);

      GenModel(const string &exstarts, const string &exends, 
            char strand, int cdsb, int cdse, const string &gi, 
            int ii, const string &si, const string &track, int trackid) throw(InvalidModel);
      ~GenModel();
      ostream& show(ostream &ous) const;
      Noschain FivePrimeUTR() const { return subchainBeforePoint(cds.begin()); }
      Noschain ThreePrimeUTR() const { return subchainAfterPoint(cds.end()); }
      int num5NoncodingExons() const { return FivePrimeUTR().numberOfRanges(); }
      int num3NoncodingExons() const { return ThreePrimeUTR().numberOfRanges(); }
      Noschain CDS() const throw (PointOutChain) { return subchain(cds); }
      Range CDSRange() const { return cds; }
      int FivePrimeUTRLength() const { return subchainBeforePoint(cds.begin()).exonLength(); }
      int ThreePrimeUTRLength() const { return subchainAfterPoint(cds.end()).exonLength(); }
      int UTRLength() const { return FivePrimeUTRLength() + ThreePrimeUTRLength(); }
      int CDSLength() const throw (PointOutChain) { return CDS().exonLength(); }
      double CDSFraction() const { return static_cast<double>(CDSLength())/length(); }

      //bool valid();
      int valid();
      bool valid(const string &gs); 
      //bool good() const;
      // information
      //string& genomicId() const { return gid; }
      const string& genomicId() const { return gid; }
      const string& modelName() const { return name; }
      const string& getName() const { return name; }
      int getId() const { return id; }
      string CDSSeq(const string &gs) const throw (PointOutChain);
      void growEnd(int icr);
      void setEnd(int ne) { Noschain::setEnd(ne); cds.setEnd(ne); }
      string toString(char sep='\t') const;
      string toSQLString(char sep=',') const;
      string toJGIModel(char sep='\t') const;
      ostream& printAllmod(ostream &ous, char sep='\t') const;
      bool operator==(const GenModel &gm) const { return Noschain::operator==(gm) && cds==gm.cds; }
      bool sameExons(const GenModel &gm) const { return Noschain::operator==(gm); }
      void makeProtein(const string &gs);
      const string& pepstr() const { return pep->getSequence(); }
      bool hasBegin() const { return hasstart; }
      bool hasEnd() const { return hasstop; }

   private:
      string gid;  
      int id;
      string name;
      Range cds;

      int cdsphase;
      int numstop;
      bool hasstart, hasstop;
      Protein *pep;
      string track;
      int trackid;

      static const int minModelLen=100;
      static const int minCDSLen=80;
};

#endif