gene.hh

/*
 * gene.hh
 *
 * License: Artistic License, see file LICENSE.TXT or 
 *          https://opensource.org/licenses/artistic-license-1.0
 */
 
#ifndef _GENE_HH
#define _GENE_HH
 
// project includes
#include "motif.hh"
#include "pp_scoring.hh"
#include "hints.hh"
 
 
// Forward declarations
 
class State;
 
class Transcript;
 
class Gene;
 
class AltGene;
 
class AnnoSequence;
 
template<class T>
struct ptr_comparison
{
    bool operator()(T* a, T* b) { return *a < *b; } 
};
 
 
class SrcEvidence {
public:
    SrcEvidence(string srcname){
    this->srcname = srcname;
    freq = 0;
    }
    string srcname;
    int freq;
    list<string> groupnames;
};
 
bool operator<(const SrcEvidence& e1, const SrcEvidence& e2);
 
class Evidence {
public:
    Evidence(bool withNames) {
    numEvidence = 0;
    this->withNames = withNames;
    }
    ~Evidence() {
    }
    void add(string source, string name);
    void add(string source);
    void print();
    list<SrcEvidence> sourceEvidence;
    int numEvidence;
    bool withNames;
};
 
Evidence *stateEvidenceSummary(State *state1, State *state2 = NULL);
 
 
class State {
public:
  long     begin;
  long     end;
  State *next;
  StateType type;
  Double prob;
  bool hasScore;
  float apostprob; // score = posterior probability 
  int sampleCount;
  Evidence *evidence;
  char truncated; // truncated left, truncated right?
  int framemod; // reading frame modifier, frame = frame given by type + framemod (for truncation)
    
  State(long first, long last, StateType t) :
      begin(first), end(last),
      next(NULL),
      type(t),
      prob(0),
      hasScore(false),
      apostprob(0),
      sampleCount(1),
      evidence(NULL),
      truncated(0),
      framemod(0)
  {}
 
  State() : 
      begin(0), 
      end(0), 
      next(NULL), 
      type(TYPE_UNKNOWN),
      prob(0),
      hasScore(false),
      apostprob(0),
      sampleCount(1),
      evidence(NULL),
      truncated(0),
      framemod(0) {}
  ~State(){
    if (evidence)
    delete evidence;
  }
  State *cloneStateSequence(){
    State *erg = new State(*this);
    if (next)
      erg->next = next->cloneStateSequence();
    else 
      erg->next = NULL;
    return erg;
  }
  // copy constructor
    State(const State& other) :
    begin(other.begin),
    end(other.end),
    next(other.next),
    type(other.type),
    prob(other.prob),
    hasScore(other.hasScore),
    apostprob(other.apostprob),
    sampleCount(other.sampleCount),
    evidence(NULL),
    truncated(other.truncated),
    framemod(other.framemod)
    {
    if (other.evidence)
        evidence = new Evidence(*other.evidence);
    }
  int frame(); // reading frame, in case it is a (coding) exon, usually depends only on the type
  bool frame_compatible(const Feature *hint);
  void addEvidence(string srcname) {if (!evidence) evidence = new Evidence(false); evidence->add(srcname);}
  bool operator< (const State &other) const;
  bool operator== (const State &other) const;
  int length() {
      return (int) (end - begin + 1);
  }
  Strand strand() {return isOnFStrand(type)? plusstrand : minusstrand;}
  State *getBiologicalState();
  void setTruncFlag(long end, long predEnd, long dnalen);
  void includeFrameModIntoType();
};
 
bool frame_compatible(State *ex1, State* ex2);
 
class StatePath {
public:
    StatePath(){
    first = NULL;
    pathemiProb = 0.0;
    intron_d = 0;
    seqname = "";
    }
    ~StatePath(){
    State *tmp;
    while( first ){
        tmp = first->next;
        delete first;
        first = tmp;
    }
    }
 
    void push(State *st){
    st->next = first;
    first = st;
    }
 
    int size() {
    int n = 0;
    for (State* st = first; st != NULL; st = st->next) {
        n++;
    }
    return n;
    }
 
    void print();
    static StatePath* condenseStatePath(StatePath *oldpath);
    Transcript* projectOntoGeneSequence(const char *genenames);
    static StatePath* getInducedStatePath(Transcript *genelist, long dnalen, bool printErrors=true);
    void reverse();
    bool operator== (const StatePath &other) const;
    bool operator< (const StatePath &other) const;
private:
    void pushIntron(long begin, long end, int frame, bool onFStrand);
 
public:
    string seqname;
    State* first;
    Double pathemiProb;
    int intron_d;
    list<PP::Match> proteinMatches;
};
 
 
int lenStateList(State *head);
 
class Transcript {
public:
    Transcript() {
    exons = introns = (State*) NULL;
    next = NULL;
    transstart = transend = -1;
    seqname = id = source = "";
    strand = plusstrand;
    complete = true;
    apostprob = 0.0;
    hasProbs = false;
    throwaway = false;
    viterbi = true;
    }
    Transcript (const Transcript& other);
    Transcript& operator= (const Transcript& other);
    virtual ~Transcript(){
    State *st, *tmp;
    list<State*> sl = getExInInHeads();
    for (list<State*>::iterator it = sl.begin(); it != sl.end(); ++it){
        st = *it;
        while( st ){
        tmp = st->next;
        delete st;
        st = tmp;
        }
    }
    }
    virtual Transcript* clone() {return new Transcript(*this);}
    Transcript *cloneGeneSequence(){
    Transcript *res = clone(); // returns a Transcript or a Gene as appropriate
    if (next)
        res->next = next->cloneGeneSequence();
    else 
        res->next = NULL;
    return res;
    }
    static void destroyGeneSequence(Transcript *head) {
    Transcript* nextHead = head;
    while (nextHead) {
        head = nextHead;
        nextHead = nextHead->next;
        delete head;
    }
    }
    void addStatePostProbs(float p);
    void setStatePostProbs(float p);
    void addSampleCount(int k);
    void setSampleCount(int k);
    virtual long geneBegin() const { return transstart;}
    virtual long geneEnd() const { return transend;}
    virtual bool isCoding() const { return false; }
    bool operator< (const Transcript &other) const;
    bool operator== (const Transcript &other) const;
    void normPostProb(float n);
    void updatePostProb(Transcript* other);
    virtual list<State*> getExInHeads() const {
    list<State*> L;
    L.push_back(exons);
    L.push_back(introns);
    return L;
    }
    virtual list<State*> getExInInHeads() const { return getExInHeads();}
    double meanStateProb();
    char* getExonicSequence(AnnoSequence *annoseq = NULL,
                bool noOffset = false) const; // CDS or whole RNA, respectively
    virtual void shiftCoordinates(long d);
    virtual bool almostIdenticalTo(Transcript *other);
    virtual void printCodingSeq(AnnoSequence *annoseq) const {}; // print nothing
    virtual void printProteinSeq(AnnoSequence *annoseq) const {};// for noncoding
    virtual void printBlockSequences(AnnoSequence *annoseq) const {}; // genes
    virtual void printGFF() const;
    virtual void printEvidence() const {}; // implemented only for coding genes
    void setStateHasScore(bool has);
    static Transcript* getGenesOnStrand(Transcript* genes, Strand strand);
    static void filterTranscriptsByMaxTracks(list<Transcript*> &gl, int maxTracks);
    virtual double supportingFraction(HintGroup *group) {return 0.0;}
public:
    State*      exons; // the exons (not UTR)
    State*      introns; // the introns between 'exons'
    long        transstart, transend; // transcription boundaries, -1 if not known
    Transcript* next;
    string  id;       // transcript id
    string  seqname;
    string  source;
    Strand  strand;
    Boolean complete; //coding region is complete
    string  geneid;   // id of the AltGene
    float   apostprob;
    bool    hasProbs;
    bool    throwaway;
    bool    viterbi;
 
    static bool gff3;
    static bool print_tss;
    static bool print_tts;
};
 
 
void filterGenePrediction(list<Transcript*> &gl, list<Transcript*> &filteredTranscripts, const char *seq, Strand strand, bool noInFrameStop, bool &hasInFrameStop, double minmeanexonintronprob=0.0, double minexonintronprob=0.0);
 
class Gene : public Transcript {
public:
    Gene() {
    utr5exons = utr3exons = utr5introns = utr3introns = (State*) 0;
    length   = clength = 0;
    weight   = 1;
    strand   = plusstrand;
    frame    = 0;
    complete5utr = complete3utr = true;
    codingstart = codingend = -1;
    supportingEvidence = incompatibleEvidence = CDSexonEvidence = CDSintronEvidence = UTR5stateEvidence = UTR3stateEvidence = NULL;
    }
    Gene(const Gene& other);
    virtual Gene* clone() { return new Gene(*this); }
    ~Gene(){ // this calls ~Transcript implicitly, which deletes exons and introns
    if (supportingEvidence)
        delete supportingEvidence;
    if (incompatibleEvidence)
        delete incompatibleEvidence;
    if (CDSintronEvidence)
        delete CDSintronEvidence;
    if (CDSexonEvidence)
        delete CDSexonEvidence;
    if (UTR5stateEvidence)
        delete UTR5stateEvidence;
    if (UTR3stateEvidence)
        delete UTR3stateEvidence;
    }
    list<State*> getExInHeads() const { list<State*> L; L.push_back(exons); L.push_back(introns); L.push_back(utr5exons); L.push_back(utr3exons); return L;}
    list<State*> getExInInHeads() const { list<State*> L = getExInHeads(); L.push_back(utr5introns); L.push_back(utr3introns); return L;}
    bool hasInFrameStop(AnnoSequence *annoseq) const;
    // void computeBC(char *seq);
    int numExons() const;
    State *lastExon() const;
    bool identicalCDS(Gene *other);
    using Transcript::almostIdenticalTo; // not really, but to prevent the compiler warning on MAC
    virtual bool almostIdenticalTo(Gene *other);
    void shiftCoordinates(long d);
    long geneBegin() const { return (transstart>=0)? transstart : codingstart;}
    long geneEnd() const { return (transend>=0)? transend : codingend;}
    virtual bool isCoding() const { return true; }
    void addUTR(State *mrnaRanges, bool complete_l=true, bool complete_r=true);
    void compileExtrinsicEvidence(list<HintGroup>  *groupList);
    double supportingFraction(HintGroup *group);
    void addSupportedStates(HintGroup *group);
    double getPercentSupported() const;
    long getCDSCoord(long loc, bool comp) const;  
    bool completeCDS() const;
    void print();
    void printGFF() const;
    void printCodingSeq(AnnoSequence *annoseq) const;
    void printProteinSeq(AnnoSequence *annoseq) const;
    void printBlockSequences(AnnoSequence *annoseq) const;
    virtual void printEvidence() const;
    void truncateMaskedUTR(AnnoSequence *annoseq);
 
    static void init();
 
    State*  utr5exons;
    State*  utr3exons;
    State*  utr5introns;
    State*  utr3introns;
    long    codingstart, codingend; // transstart <= codingstart <= codingend <= transend, if not -1
    bool    complete5utr;
    bool    complete3utr;
  
    int     length;
    int     clength;
    int     frame;
    BaseCount bc;
    int weight;
    Evidence *supportingEvidence;
    Evidence *incompatibleEvidence;
    Evidence *CDSexonEvidence;
    Evidence *CDSintronEvidence;
    Evidence *UTR5stateEvidence;
    Evidence *UTR3stateEvidence;
 
    list<PP::Match> proteinMatches;  // true if gene matches protein profile
 
    static bool print_start;
    static bool print_stop;
    static bool print_introns;
    static bool print_cds;
    static bool print_exonnames;
    static bool stopCodonExcludedFromCDS;
    static bool print_utr;
    static bool print_blocks;
};
 
 
class AltGene {
public:
    list<Transcript*> transcripts;
    long mincodstart; // leftmost position of any start codon (if coding)
    long maxcodend; //
    Strand strand;
    string id;
    string seqname;
    float apostprob;
    bool hasProbs;
 
    AltGene() {
    mincodstart = maxcodend = -1;
    strand = plusstrand;
    apostprob = 0.0;
    hasProbs = 0;
    }
    bool operator< (const AltGene &other) const;
    void addGene(Transcript* tx);
    bool overlaps(Transcript *tx);
    void shiftCoordinates(long d);
    void sortTranscripts(int numkeep=-1);
    void deleteSuboptimalTranscripts(bool uniqueCDS);
    long minTransBegin();
    long  maxTransEnd();
    bool isCoding(){ return transcripts.empty() || dynamic_cast<Gene*> (transcripts.front());}
};
 
Transcript* getPtr(list<AltGene> *gl);
 
void printGeneList(list<AltGene> *genelist, AnnoSequence *annoseq, bool withCS, bool withAA, bool withEvidence);
void printGeneList(Transcript* seq, AnnoSequence *annoseq, bool withCS, bool withAA);
void printGeneSequence(Transcript* seq, AnnoSequence *annoseq = NULL, bool withCS=false, bool withAA=true);
list<Gene*>* sortGenePtrList(list<Gene*>);
list<AltGene> *reverseGeneList(list<AltGene> *altGeneList, long endpos);
list<AltGene>* groupTranscriptsToGenes(list<Transcript*> &transcripts);
 
void reverseGeneSequence(Transcript* &seq, long endpos);
void postProcessGenes(list<AltGene> *genes, AnnoSequence *annoseq);
Gene* promoteToCoding(Transcript* tx, AnnoSequence *as); // make a coding gene from a non-coding if it contains a good ORF
 
class Annotation {
public:
  Annotation() {
    genes = lastGene = (Gene*) 0;
    forwardGenes =  backwardGenes  = lastForwardGene = lastBackwardGene = (Gene*) 0;
    forwardPath = condensedForwardPath = backwardPath = condensedBackwardPath = (StatePath*) 0;
    path = condensedPath = (StatePath*) 0;
    emiProb = 1.0;
    forwardEmiProb = 1.0;
    backwardEmiProb = 1.0;
  }
  ~Annotation() {
    Transcript::destroyGeneSequence(genes);
    if (path)
      delete path;
    if (condensedPath)
      delete condensedPath;
 
    Transcript::destroyGeneSequence(forwardGenes);
    Transcript::destroyGeneSequence(backwardGenes);
    if (forwardPath)
      delete forwardPath;
    if (backwardPath)
      delete backwardPath;
    if (condensedForwardPath)
      delete condensedForwardPath;
    if (condensedBackwardPath)
      delete condensedBackwardPath;
  }
 
  void appendGene(Transcript* gene);
  void appendForwardGene(Transcript* gene);
  void appendBackwardGene(Transcript* gene);
  const void printGFF() const;
 
  StatePath *path;
  StatePath *condensedPath;
 
  StatePath *forwardPath;            
  StatePath *backwardPath;           
  StatePath *condensedForwardPath;   
  StatePath *condensedBackwardPath;  
 
  Transcript *genes;
  Transcript *forwardGenes;  
  Transcript *backwardGenes; 
 
  Double emiProb;
  Double forwardEmiProb;  
  Double backwardEmiProb; 
 
private:
  Transcript *lastGene;
  Transcript *lastForwardGene, *lastBackwardGene; 
};
 
 
class AnnoSequence {
public:
    AnnoSequence(){
    length = 0;
    sequence = 0;
    seqname = 0;
    next = (AnnoSequence*) 0;
    anno = (Annotation*) 0;
    weight = 1;
    offset = 0;
    }
    ~AnnoSequence(){
    if (sequence)
        delete[] sequence;
    if (seqname)
        delete [] seqname;
    if (anno)
        delete anno;
    }
    static void deleteSequence(AnnoSequence *head){
    AnnoSequence* nextHead = head;
    while(nextHead) {
        head = nextHead;
        nextHead = nextHead->next;
        delete head;
    }
    }
    void setWeight(int w) {
    weight = w;
    if (!anno)
        return;
    
    for (Transcript* t = anno->genes; t != NULL; t = t->next){
        Gene *g = dynamic_cast<Gene*> (t);
        if (g)
        g->weight = w;
    }
    }
    /*
     * reverse complement the sequence and all genes
     */
    AnnoSequence* getReverseComplement();
    void printGFF();
 
    char*        seqname;
    long         length;
    long         offset; // if >0 this is the number of bases that have been cut out in the front of the original input sequence
    char*        sequence;
    BaseCount    bc;
    AnnoSequence *next;
    Annotation*  anno;
    int          weight;
};
 
bool isAllLC(AnnoSequence *seqs);
 
 
/*
 * Print a warning if --softmasking=1 (default) and all sequences are
 * completely in lowercase characters. This can happen in particular
 * if the input is from GenBank.
 */
void warnAllLowerCase(AnnoSequence *annoseq);
 
 
 
class AnnoSeqGeneIterator {
public:
  AnnoSeqGeneIterator(const AnnoSequence* annoseqHead){
    annoseq = annoseqHead;
    while (annoseq && !(annoseq->anno && annoseq->anno->genes))
      annoseq = annoseq->next;
    if (annoseq)
      gene = annoseq->anno->genes;
    else 
      gene = NULL;
  };
  friend const AnnoSeqGeneIterator& operator++(AnnoSeqGeneIterator& gi);
    
  bool hasMoreElements(){
    return (annoseq != NULL && gene != NULL);
  }
  const AnnoSequence *annoseq;
  const Transcript *gene;
};
 
class CompareStatePathPtr{
public: 
  CompareStatePathPtr(){ }
  ~CompareStatePathPtr(){ }
  bool operator()(StatePath *first, StatePath *second){
    return (first->pathemiProb > second->pathemiProb);
  }
};
 
class StatePathCollection {
public:
  StatePathCollection(){
    sorted = false;
  };
  ~StatePathCollection(){
    StatePath *p;
    while (!pathlist.empty()){
      p = pathlist.front();
      pathlist.pop_front();
      delete p;
    }
  }
  void addPath(StatePath *p){
    if (!containsPath(p))
      pathlist.push_back(p);
    sorted = false;
  }
  int size(){
    return pathlist.size();
  }
  void sort(){
    pathlist.sort(CompareStatePathPtr());
    sorted = true;
  }
  void printAPosterioriProbs(ostream& out) {
    for(list<StatePath*>::iterator it=pathlist.begin(); it!=pathlist.end(); it++)
      out << (*it)->pathemiProb << " , ";
    out << endl;
  }
  bool containsPath(StatePath *p);
  int positionInCollection(StatePath *p);
 
private:
  list<StatePath*> pathlist;
  bool sorted;
};
 
void examineBaseCountOfGeneSeq(AnnoSequence  *as);
 
 
class FreqSegment{
public:
    FreqSegment(long s, int f){
    start = s;
    freq = f;
    }
    long start;
    int freq;
};
 
#endif   //  _GENE_HH