RNAModel.h

Go to the documentation of this file.

#ifndef RNAMODEL_H
#define RNAMODEL_H

#include "GenModel.h"
#include <string>
#include "bioseq.h"
#include "Coverdepth.h"
#include "strformat.h"
#include <cassert>
#include <cmath>

using namespace std;

class OutsideGenomicSequence : public exception {
   private:
      string message;
   public:
      OutsideGenomicSequence(const string& msg) throw() 
         : message(msg) { }
      OutsideGenomicSequence() throw() : message("Position outside Genomic Length") { }
      ~OutsideGenomicSequence() throw() { }
      const char* what() const throw() {
         return message.c_str();
      }
};


enum modeltype { RNA, mRNA, EST, JGI};

/*
class Intronbound {
   public:
      Intronbound() : seq(), count(0) { }
      Intronbound(const char ss[5], int cc) : seq(ss), count(cc) { }
      Intronbound(const Intronbound& ib) : seq(ib.seq), count(ib.count) { }

      char seq[5];
      int count;
};
*/

class RNAModel : public Noschain {
   public:
      RNAModel() : Noschain(), genid(0) { oid=nextId(); }

      explicit RNAModel(const Noschain &seg, const string &chrom, const string &gs)
         : Noschain(seg), gid(chrom), gseq(&gs), genid(0), oid(0) { 
            oid=nextId(); subsequence(rna, *gseq);
      }
      RNAModel(const string &exstr, const string &gname, const string &gs, 
            int oi)
         : Noschain(exstr), gid(gname), gseq(&gs), oid(oi), genid(0) {
            subsequence(rna, *gseq); }
      RNAModel(const string &exstr, const string &gname, const string &gs, 
            int oi, int gene_id)
         : Noschain(exstr), gid(gname), gseq(&gs), oid(oi), genid(gene_id) {
            subsequence(rna, *gseq); }
      RNAModel(const string &exstarts, const string &exends,
            char strand, const string& gi, const string& genomic)
         throw(Badinput, exception)
         : Noschain(exstarts, exends, strand), gid(gi), gseq(&genomic),
            genid(0), oid(serialid++) {
               if (exons[0].isNull() || exons.back().isNull()) 
                  throw Badinput("Null end exon");
            subsequence(rna, *gseq);
      }
      RNAModel(const string &exstarts, const string &exends,
            char strand, const string& gi, const string& genomic,
            int oo)
         throw(Badinput)
         : Noschain(exstarts, exends, strand), gid(gi), gseq(&genomic),
            genid(0), oid(oo) {
               if (exons[0].isNull() || exons.back().isNull()) 
                  throw Badinput("Null end exon");
            subsequence(rna, *gseq);
         }
      RNAModel(const RNAModel &rm) : Noschain(rm), gid(rm.gid),
         gseq(rm.gseq), rna(rm.rna), oid(serialid++), genid(rm.genid) { }
      RNAModel& operator=(const RNAModel &rm);
      void assignNewId() { oid=nextId(); }
      int getOid() const { return oid; }
      void setOid(int oi) { oid=oi; }
      void setNextOid() { oid=nextId(); }
      void assignNewGeneId() { genid=nextGeneId(); }
      void modelidAsGeneid() { genid=oid; }
      void assignCurrentGeneId() { genid=currentGeneId(); }
      int getGeneId() const { return genid; }
      int geneId() const { return genid; }
      void setGeneId(const int gene_id) { genid=gene_id; }
      void setGenomic(const string &gs) { gseq=&gs; }
      void setGenomicPointer(const string *gs) { gseq=gs; }
      void trimAfterPoint(const int p) throw (PointOutChain);
      void trimBeforePoint(const int p) throw (PointOutChain);
      const string& RNASequence() const { return rna; }
      int genomicIndex(int mrnaidx) const throw (OutsideGenomicSequence);
      int RNAIndex(int gidx) const throw (PointOutChain);
      bool sameExons(const RNAModel &gm) const { return Noschain::operator==(gm); }
      const string* getGenomicPointer() const { return gseq; }
      int exonLength() const { return rna.length(); }
      virtual bool valid() const;
      //bool operator==(const RNAModel &rmd) const { return rna == rmd.rna; }
      bool sameRNA(const RNAModel &rmd) const { return rna==rmd.rna; }
      const string& genomicId() const { return gid; }
      int RNALength() const { return rna.length(); }
      const string& RNAString() const { return rna; }
      void resetRNA() throw(Badinput) { subsequence(rna,*gseq); }
      void setRNA(const string &rseq) { rna=rseq; }
      virtual void reset() { subsequence(rna,*gseq); }
      void setGenomicId(const string &gg) { gid=gg; }
      bool append(Noschain &mod);
      RNAModel& reverse();
      
      string seqGenomic() const;
      int genomicLength() const { return gseq->length(); }
      virtual modeltype objtype() const { return RNA; }
      ostream& writeExon(ostream &ous, char sep='\t') const;
      ostream& print(ostream& ous) const;
      void intronBound(map<string,int> &spsite) const;

      static const char exonheader[];

      static int currentId() { return serialid; }
      static int nextId() { return serialid++; }
      static void advanceId() { ++serialid; }
      static void rollbackId() { --serialid; }
      static void setId(int oi) { serialid=oi; }

      static int nextGeneId() { return geneid++; }
      static int currentGeneId() { return geneid; }
      static void advanceGeneId() { ++geneid; }

   protected:

      string gid;
      const string *gseq;
      string rna;
      int oid; 
      int genid; 

      static int serialid; 
      static int exonid;
      static int geneid;
};


class mRNAModel : public RNAModel {
   public:
      mRNAModel() : RNAModel(), cdsb(-1), cdse(-1),
         gcdsb(-1), gcdse(-1), pep(), frame(0) { }
      explicit mRNAModel(const Noschain &seg, const string &chrom, const string &gs);
      explicit mRNAModel(const Noschain &seg, int gcdsb_, int gcdse_, const string &chrom, const string &gs, int fr=0);
      mRNAModel(const string &exstr, const string &gname, const string &gs,
            int oi, int gcb, int gce, int fr=0);
      mRNAModel(const string &exstr, const string &gname, const string &gs,
            int oi, int gcb, int gce, int gene_id, int fr=0);
      mRNAModel(const string &exstarts, const string &exends, char strand, 
            int gcb, int gce, const string& gi, const string& genomic)
         throw(PointOutChain, Badinput, exception);
      mRNAModel(const string &exstarts, const string &exends, char strand, 
            int gcb, int gce, const string& gi, const string& genomic, int oo)
         throw(PointOutChain, Badinput);

      mRNAModel(const mRNAModel &mm) : RNAModel(mm), cdsb(mm.cdsb),
         cdse(mm.cdse), gcdsb(mm.gcdsb), gcdse(mm.gcdse), 
         pep(mm.pep), frame(mm.frame) { }
      char guessStrand() const;

      mRNAModel& operator=(const mRNAModel &mm);
      bool samePeptide(const mRNAModel &mm) const { return pep==mm.pep; }
      bool sameGene(const mRNAModel &mod) const;
      Noschain FivePrimeUTR() const { return subchainBeforePoint(gcdsb); }
      Noschain ThreePrimeUTR() const { return subchainAfterPoint(gcdse); }
      int num5NoncodingExons() const { return FivePrimeUTR().numberOfRanges(); }
      int num3NoncodingExons() const { return ThreePrimeUTR().numberOfRanges(); }
      Noschain CDSChain() const throw (PointOutChain); //{ return subchain(gcdsb, gcdse); }
      Range CDSRange() const { return Range(gcdsb, gcdse); }
      pair<int,int> genomicCDSBound() const { return pair<int,int>(gcdsb, gcdse); }
      Range RNACDSRange() const { return Range(cdsb,cdse); }
      pair<int,int> RNACDSBound() const { return pair<int,int>(cdsb, cdse); }
      int CDSLength() const { return cdse-cdsb+1; }
      int genomicCDSLength() const { return abs(gcdse-gcdsb)+1; }
      int genomicCDSEnd() const { return gcdse; }
      int genomicCDSBegin() const { return gcdsb; }
      const string& proteinSequence() const { return pep; }
      int proteinLength() const { return pep.length(); }
      int proteinLengthNoTail() const;
      int getFrame() const { return frame; }
      string CDSSeq() const { return rna.substr(cdsb-1, cdse-cdsb+1); }
      string CDSSequence() const;
      char CDSDirection() const { if (gcdsb < gcdse) return '+'; return '-'; }
      double CDSFraction() const { return static_cast<double>(CDSLength())/length(); }
      double CDSFractionRNA() const { return static_cast<double>(CDSLength())/exonLength(); }
      double CDSFractionGenomic() const { return static_cast<double>(CDSLength())/length(); }
      int FivePrimeUTRLength() const { return cdsb-1; }
      int ThreePrimeUTRLength() const { return exonLength() - cdse; }
      int UTRLength() const { return FivePrimeUTRLength() + ThreePrimeUTRLength(); }
      void UTR3Sequence(string &seq) const { if (cdse==rna.length()) seq=""; else seq=rna.substr(cdse); }
      string UTR3Sequence() const { if (cdse==rna.length()) return ""; return rna.substr(cdse); }
      bool hasStart() const { return pep[0] == 'M'; }
      bool hasStop() const { return pep[pep.length()-1] == '*'; }
      bool complete() const { return hasStart() && hasStop(); }
      //bool operator==(const mRNAModel &md) const { 
       //  return Noschain::operator==(md) && pep==md.pep; }
      int numberOfInternalStops() const;
      const string& getProtein() const { return pep; }
      void trimAfterPoint(const int p) throw (PointOutChain);
      void trimBeforePoint(const int p) throw (PointOutChain);
      void trimBeforePoint(const int gp, const int rp) throw (PointOutChain);
      void setProtein(const string &pseq) { pep=pseq; }
      void setRNACDS(int bb, int ee) throw(OutsideGenomicSequence);
      void resetRNACDS(int bb, int ee) throw(OutsideGenomicSequence) {
         setRNACDS(bb,ee); resetProtein(); }
      void setGenomicCDS(int gb, int ge) throw(OutsideGenomicSequence, InvalidModel);
      void setGenomicCDS(pair<int, int> gcr) {
        //try { setGenomicCDS(gcr.first, gcr.second); } catch (exception &err) { throw err; } }
      setGenomicCDS(gcr.first, gcr.second); assert(cdsb<cdse); }
      bool growCDS3Prime(int len) throw(OutsideGenomicSequence);
      bool trimCDSTail();
      bool trimCDSStop();
      void reset();
      void resetProtein() { translate(pep, rna, cdsb+frame, cdse); }
      void setLongestCDSAndProtein();
      bool append(mRNAModel &mod, int &comment);
      string JGIModelRow(const char sep='\t') const;
      ostream& printJGIModelRow(ostream &ous, const char sep='\t') const;
      string JGITranscriptRow(const char sep='\t') const;
      ostream& printJGITranscriptRow(ostream &ous, const char sep='\t') const;
      ostream& printJGITranscriptRowNoId(ostream &ous, const char sep='\t') const;
      string JGIProteinRow(char sep='\t') const;
      //string JGILinkToAssRow(char sep='\t') const;
      string sfCDSGenomic() const;
      Noschain CDSOnewayChain() const;
      void CDSOnewayChain(Noschain& ch) const;
      string sfCDSTranscript() const {
         return (CDSChain().onewayChain().exonOnlyChain()+(cdsb-1)).asDelimitedString(':', ','); 
      }
      string sfExonsProtein() const;
      ostream& show(ostream &ous) const;
      bool genuine() const
         { return pep.length()>50 && complete()
               && ThreePrimeUTRLength()>30 && ThreePrimeUTRLength() < 1800 
               && FivePrimeUTRLength()>30 && FivePrimeUTRLength()< 1800 
               && (CDSFractionRNA() > 0.5 || CDSLength() > 1000) 
               && num3NoncodingExons() < 3 && num5NoncodingExons() < 3; }
      bool semiGenuine() const { return num3NoncodingExons() < 3 && num5NoncodingExons() < 3 && CDSLength() > 140 && complete() && ThreePrimeUTRLength() < utrlen3max && FivePrimeUTRLength() < utrlen5max; }
      bool isStar() const { return (CDSFraction() > 0.6 || (ThreePrimeUTRLength() < utrlen3max && FivePrimeUTRLength() < utrlen5max && CDSLength() > 240)) && num3NoncodingExons()<3 && num5NoncodingExons() < 3 && complete(); }
      modeltype objtype() const { return mRNA; }
      bool valid() const;
      ostream& print(ostream &ous) const;
      mRNAModel& reverse(int newRNACDSB, int newRNACDSE);
      string name(const string &prefix) { return prefix + itos(getOid()); }

      static void setShortestModel(int len) { shortestmodel=len; }
      void writetab(ostream& mod, ostream &ex, ostream &track,
            ostream &orna, ostream &oprt, char sep='\t') const;
      ostream& writeModelTable(ostream &ous, char sep='\t') const;

      static const char modelheader[];
      static const char jgiModelCol[];
      static const char jgiTranscriptCol[];
      static const char jgiProteinCol[];

      static int shortestpep;
      static int shortestmodel;
      static int utrlen5max, utrlen3max;
      
   protected:
      int cdsb, cdse;
      int gcdsb, gcdse; 
      string pep;
      int frame; // 0, 1, 2

       static char header[350];
};


class ESTAssembly;

class mRNAModelUpdate : public mRNAModel {
   public: 
      explicit mRNAModelUpdate(const mRNAModel &mm, int ec) 
         : mRNAModel(mm), numest(ec), comment(0) { }
      mRNAModelUpdate(const mRNAModel &mm) 
         : mRNAModel(mm), numest(0), comment(0) { }
      mRNAModelUpdate(const mRNAModelUpdate &um) : mRNAModel(um), numest(um.numest), estcover(um.estcover), comment(um.comment) { }
      void addESTCover(const mRNAModel &mm, const char direction);
      void addESTCover(const ESTAssembly &mm);
      void setESTCount(int ec) { numest=ec; }
      void addESTCount(int ec) { numest += ec; }
      void writetab(ostream& mod, ostream &ex, ostream &track,
            ostream &orna, ostream &oprt, char sep='\t') const;
      ostream& writeModelTable(ostream &ous, char sep='\t') const;

      string name() const { return "MIX" + itos(getOid()) + "_" + itos(numest) + "_" + itos(static_cast<int>(floor(ESTCoverage()*100+0.5))); }
      string title() const { return itos(numest) + " EST updated Ab initio model"; }
      bool append(ESTAssembly &mod);
      float ESTCoverage() const { return static_cast<float>(estcover.exonLength())/exonLength(); }
      ostream& show(ostream &ous) const;

      static const char modelheader[];
      static const char* jgiModelColumns();
      static const char* jgiTranscriptColumns();
      static const char* jgiProteinColumns();


   private:
      int numest;
      RangeChain estcover;
      int comment; 
};

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

class JGIModel : public mRNAModel {
   public:
      JGIModel(const string &exstarts, const string &exends,
            char strand, int cdsb, int cdse, const string &gi,
            int ii, const string &si, const string &genomic) throw(Badinput, InvalidJGIModel, InvalidModel, PointOutChain, exception);
      JGIModel(const string &exstarts, const string &exends,
            char strand, int cdsb, int cdse, const string &gi,
            const string &si, const string &genomic) throw(Badinput, InvalidJGIModel, InvalidModel, PointOutChain, exception);
      bool valid();
      string toJGIString(char sep='\t') const;
      string toSQLString(char sep=',') const;
      string toString(char sep='\t') const;
      void growEnd(int icr);
      int getId() const { return getOid(); }
      const string& getName() const { return name; }
      ostream& show(ostream &ous) const;

      ostream& print(ostream &ous) const;
      modeltype objtype() const { return JGI; }

   private:

      /* RNAModel class has oid. This id should be
       * removed. Eliminating this member, use oid from
       * RNAModel.
       */
      //int id; 
      string name;
      // using frame from mRNA parent class
      // removing cdsphase
      //int cdsphase; // is this the same as frame from mRNAModel?
      //int numstop;
};

bool useFullORF(int lenfull, int lenpartial);

class ESTAssembly : public mRNAModel {
   public:
      ESTAssembly() : mRNAModel(), congid(-1), numest(0), basecovdep(0), comment(), intronfix(0), relprofh(1) { }
      explicit ESTAssembly(const Noschain &seg, int nume, Coverdepth *bc,
            const string &chrom, const string &gs, int cngid);

      explicit ESTAssembly(const Noschain &seg, int gcdsb_, int gcdse_, int nume, 
            Coverdepth *bc, const string &chrom, const string &gs, int cngid, int fr);
      explicit ESTAssembly(const Noschain &seg, int nume, 
            const string &chrom, const string &gs, int cngid);
      ESTAssembly(const string &exons, const string &genomicid,
            const string &genomicseq, int modelid, int gcdsB, int gcdsE, 
            int congregationid, int estcount, int gnid, int phase, float prof2star)
         : mRNAModel(exons, genomicid, genomicseq, modelid, 
               gcdsB, gcdsE, gnid, phase),
            congid(congregationid), numest(estcount), basecovdep(0),
            comment(), intronfix(0), relprofh(prof2star)
            { setNameAndTitle(); }
      ESTAssembly(const ESTAssembly &ea) : mRNAModel(ea), 
         congid(ea.congid), numest(ea.numest), estname(ea.estname),
         esttitle(ea.esttitle), pepname(ea.pepname), peptitle(ea.peptitle), comment(ea.comment), intronfix(ea.intronfix), basecovdep(ea.basecovdep), relprofh(ea.relprofh) { }
      ~ESTAssembly() { if (basecovdep != 0) delete basecovdep; }
      void setNameAndTitle();
      void setName();
      void setTitle();
      void setTitle(const string &prefix);
      void addTitleTag(const string &tag);
      ESTAssembly& operator=(const ESTAssembly &ea);
      const string& name() const { return estname; }
      void write(ostream& mod, ostream &exon, ostream &track,
            ostream &ousest, ostream &ouspep, char sep='\t') const;
       void writetab(ostream& mod, ostream &exon, ostream &track,
            ostream &ousest, ostream &ouspep, char sep='\t') const;
       void writeModel(ostream &ous, char sep='\t') const;
       float RNACodingFraction() { return CDSLength()/(float)RNALength(); }
       //bool breakSuffixModel(ESTAssembly &mod);
       ESTAssembly* breakSuffixModel() throw(PointOutChain);
       ESTAssembly* breakPrefixModel() throw(PointOutChain);
       ESTAssembly* prune3PrimeUTR() throw (PointOutChain);
       ESTAssembly* prune5PrimeUTR() throw (PointOutChain);
      /* break up current model into smaller models.
        * If not successful then the vector will be
        * empty. Try to find optimal ORFs with the 
        * RNA, then produce new models.
        * This method is independent of the break-prefix, suffix methods
        * They share very little, and used different algorithms.
      vector<ESTAssembly*> breakup();
      I have moved it to ESTAssemblyid class since this
      one was never used directly.
        */
      bool shouldBreakSuffix() const;
      bool shouldBreakPrefix() const;
      bool isChimera() const;
      modeltype objtype() const { return EST; }
      ostream& show(ostream &ous) const;
      ostream& print(ostream &ous) const;
      void showBaseCoverInfo(ostream &ous=cerr) const { 
         if (basecovdep != 0) ous << *basecovdep << endl; }
      int ESTCount() const { return numest; }
      int numberOfEST() const { return numest; }
      bool validProfile() const;
      void addComment(const string &comm);
      void writeCommentTab(ostream &ous) const { 
         if (!comment.empty()) ous << getOid() << '\t' << comment << endl; }
      bool hasComment() const { return !comment.empty(); }
      bool noComment() const { return comment.empty(); }
      void fixIntronBound();
      int intronFixState() const { return intronfix; }
      void checkIntronBound();
      bool genuine() const {
         return mRNAModel::genuine() && intronfix < 2; }
      bool semiGenuine() const {
         return mRNAModel::semiGenuine() && intronfix<2; }
      bool isStar() const {
         return mRNAModel::isStar() && intronfix<2; }
      bool findDip(int b, int e, Dip &dip) const;
      ostream& showProfile(ostream &ous) const { ous << *basecovdep << endl; return ous; }
      float averageProfileHeight() const { if (basecovdep == 0) return 1; else return basecovdep->averageHeight(); }
      int maxProfileHeight() const { if (basecovdep == 0) return 1; else return basecovdep->maxHeight(); }
      pair<float,float> relativeProfileHeight(const ESTAssembly &ea) const { return basecovdep->relativeHeight(*(ea.basecovdep)); }
      void assignRPH(const ESTAssembly &star) { relprofh=star.relativeProfileHeight(*this).first; }
      float getRelprofh() const { return relprofh; }
      void setRelprofh(const float rh) { relprofh=rh; }
      void releaseProfile() { basecovdep=0; }
      ESTAssembly& reverse(int RNACDSb, int RNACDSe) { mRNAModel::reverse(RNACDSb, RNACDSe); 
         if (basecovdep != 0) basecovdep->reverseStrand(); }

       static const char* JGIModelColumns();
       static const char* JGITranscriptColumns();
       static const char* JGIProteinColumns();
       static bool ORFLongEnough(pair<int,int> &bound, int ORFType, int trlen); 

       static const char modelheader[];
       //static const char exonheader[]; // move to RNA class
       static const char trackheader[];

   protected:

       int congid;
       int numest;
       Coverdepth *basecovdep; 
       //pair<float,float> relprofh;
       float relprofh;
       string estname, esttitle;
       string pepname, peptitle;
       string comment;
       int intronfix;

       ESTAssembly* budPlusSuffixModel(pair<int,int> gcbts, int rcdsb, int rcdse, int suffix_gcdsb, int suffix_gcdse, const string &pep, int rcut1, int rcut2) throw(PointOutChain);
       ESTAssembly* budPlusSuffixModel(pair<int,int> gcuts, int rcdsb, int rcdse, int suffix_gcdsb, int suffix_gcdse, const string &pep) throw (PointOutChain) {
          try {
          return budPlusSuffixModel(gcuts, rcdsb, rcdse, suffix_gcdsb, suffix_gcdse, pep, RNAIndex(gcuts.first), RNAIndex(gcuts.second)); } 
          catch (PointOutChain &pocerr) { cerr << "budPlusSuffixModel() short version problem\n"; throw pocerr; } }
       ESTAssembly* budMinusSuffixModel(pair<int,int> gcuts, int suffix_gcdsb, int suffix_gcdse, const string &subRNA, const string &pep, int usefull, pair<int,int> subcds, int rcut1, int rcut2) throw (PointOutChain);
       ESTAssembly* budMinusSuffixModel(pair<int,int> gcuts, int suffix_gcdsb, int suffix_gcdse, const string &subRNA, const string &pep, int usefull, pair<int,int> subcds) throw (PointOutChain) {
          try {
          return budMinusSuffixModel(gcuts, suffix_gcdsb, suffix_gcdse, subRNA, pep, usefull, subcds, RNAIndex(gcuts.first), RNAIndex(gcuts.second));
          } catch (PointOutChain &pocerr) {
             cerr << "budMinusSuffixModel() short version problem\n";
             throw pocerr;
          }
       }

       ESTAssembly* budPlusPrefixModel(pair<int,int> gcuts, int rcut1, int rcut2,
             int prefix_gcdsb, int prefix_gcdse, pair<int,int> subcds,
             int usefull, const string &maxPep) throw(PointOutChain);
       ESTAssembly* budMinusPrefixModel(pair<int,int> gcuts, int rcut1, int rcut2, 
             int prefix_gcdsb, int prefix_gcdse, pair<int,int> subcds, 
             const string &pep, const string &subRNA) throw(PointOutChain);

       ESTAssembly* budTemplate() const;
       void resetNumest(ESTAssembly *mod);
       void resetSuffixProfile(ESTAssembly *mod, pair<int,int> gcuts);
       void resetPrefixProfile(ESTAssembly *mod, pair<int,int> gcuts);
       void setCDSInfo();

       /* used for construction move this to base class 
        * mRNAModel*/
       //static char header[250];

       //static int id;
       //static int exonid;
       static int chimera_numestcut; // set to 100
       static int chimera_peplencut;    // set to 100
       static int chimera_partial_peplencut; // set 30
       static int orfspace; // now set to be short
};

//template<class T>
class ESTAssemblyid : public ESTAssembly {
   public:
      explicit ESTAssemblyid(const Alnchainid &seg, const string &chrom, 
            const string &gs, int cngid) 
         : ESTAssembly(seg,seg.numEST(),chrom,gs,cngid), estids(seg.getESTIds()) { }
      explicit ESTAssemblyid(const Noschain &seg, int nume, Coverdepth *bc,
            const string &chrom, const string &gs, int cngid)
         : ESTAssembly(seg,nume,bc,chrom,gs,cngid), estids() { }

      explicit ESTAssemblyid(ESTAssembly &est) : ESTAssembly(est), estids() { est.releaseProfile(); }
      explicit ESTAssemblyid(const Noschain &seg, int gcdsb_, int gcdse_, int nume, 
            Coverdepth *bc, const string &chrom, const string &gs, int cngid, int fr) : ESTAssembly(seg,gcdsb_,gcdse_,nume,bc, chrom, gs, cngid,fr), estids() { }
      //ESTAssemblyid(ESTAssembly *est) : estmod(est), estids() { }

      void assignESTId(const vector<string> &ids) { estids=ids; }

      void addESTId(const string &eid) { estids.push_back(eid); }
      void addESTId(const vector<string> &eids) { estids.insert(estids.end(), eids.begin(), eids.end()); }
      vector<string>& getESTIds() { return estids; }
      string ESTIdString(const char sep=',') const;
      ostream& outputESTIds(ostream &ous, const char* sep=",") const;

      ESTAssemblyid* breakSuffixModel() throw (PointOutChain);
      ESTAssemblyid* breakPrefixModel() throw (PointOutChain);
      vector<ESTAssemblyid*> breakup();
      void writeModel(ostream &ous, char sep='\t') const;
      void writetab(ostream& mod, ostream &exon, ostream &track,
            ostream &ousest, ostream &ouspep, char sep='\t') const;
      //int ESTCount() const { return estids.size(); }
      //int numberOfEST() const { return estids.size(); }
      // should use parent class header!
      //static const char modelheader[];
      static const char* getModelheader();

   private:

      vector<string> estids;
};

void clumpModelIntoGene(set<ESTAssemblyid*, lessByChainDirectionPtr> &mod,
        list<list<ESTAssemblyid*> > &loci);

#endif