namgene.hh

/*
 * namgene.hh
 *
 * License: Artistic License, see file LICENSE.TXT or 
 *          https://opensource.org/licenses/artistic-license-1.0
 */
 
#ifndef _NAM_GENE_HH
#define _NAM_GENE_HH
#include <climits>
 
// project includes
#include "extrinsicinfo.hh"  // SequenceFeatureCollection
#include "vitmatrix.hh"
#include "matrix.hh"
#include "pp_scoring.hh"
#include "statemodel.hh"
 
class NAMGeneError : public ProjectError {
public:
    NAMGeneError( string msg ) : ProjectError( msg ) { }
};
 
class NAMGene {
public:
    NAMGene();
    ~NAMGene() {}
 
    void readModelProbabilities( int number = 1);
    StatePath* getSampledPath(const char *dna, const char* seqname = NULL);
    StatePath* getViterbiPath(const char *dna, const char* seqname = NULL);
    StatePath* getTrainViterbiPath(const char *dna, SequenceFeatureCollection *sfc);
    Transcript* doViterbiPiecewise(SequenceFeatureCollection& sfc, AnnoSequence *annoseq, Strand strand, unsigned int& num_pieces);
 
    Double getEmissionProbability();
 
    const ViterbiMatrixType& getForwardVariables( ) {
    return forward;
    } 
    
    const ViterbiMatrixType& getViterbiVariables( ) {
    return viterbi;
    }
 
    /*
     * get the type corresponding to a state number
     * @param   i the thumber of a state, as in the properties
     * @return  the enumeration type
     */
    StateType getStateType(int i);
    
    int getStateIndex(StateType type);
    Double getPathEmiProb(StatePath *path, const char *dna, SequenceFeatureCollection& sfc, int countStart=-1, int countEnd=-1);
    void setNeedForwardTable(bool b){needForwardTable = b;}
    // set the path and emiProbs for all annotations in annoseq (for option scoreTx)
    void setPathAndProb(AnnoSequence *annoseq, FeatureCollection &extrinsicFeatures);
    void setAllTranscripts(list<Transcript*> *tl) {sampledTxs = tl;}
    list<Transcript*> *getAllTranscripts() {return sampledTxs;}
    void getPrepareModels(const char *dna, int len) {prepareModels(dna, len);}
private:
    void viterbiAndForward(const char* dna, bool useProfile=false);
    
    /*
     * repeatedly do the viterbi algorithm on pieces of dna, no strands
     */
    list<AltGene>* getStepGenes(AnnoSequence *annoseq, SequenceFeatureCollection& sfc, Strand strand, bool onlyViterbi=true);
 
    list<AltGene>* findGenes(const char *dna, Strand strand, bool onlyViterbi=true);
    int getNextCutEndPoint(const char *dna, int beginPos, int maxstep, SequenceFeatureCollection& sfc);
    void readTransAndInitProbs( );
    void readOvlpLenDist( );
    void checkProbsConsistency( );
    void computeReachableStates( );
    void createStateModels( );
    void setStatesInitialProbs();
    void initAlgorithms(); // called for new DNA
    void updateToLocalGCEach(int idx, int from = -1, int to = -1);
    void printDPMatrix(); // print Viterbi matrix
    void prepareModels(const char*dna, int len);
 
    long tryFindCutEndPoint(StatePath *condensedExamPath, int examIntervalStart, int examIntervalEnd, list<Feature> *groupGaps, bool onlyInternalIR);
private:
    ViterbiMatrixType      forward;
    ViterbiMatrixType      viterbi; 
    Matrix<Double>      transitions;
    vector<Double>      initProbs;
    vector<Double>      termProbs;
    // Whether the state is reachable at all with prob > 0
    vector<Boolean> stateReachable;
    vector<StateModel*> states;  
    Integer             statecount;
    vector<StateType>   stateMap;
    int                 lastParIndex;
    Boolean proteinOutput;
    Boolean codSeqOutput;
    Boolean noInFrameStop;
    double minmeanexonintronprob;
    double minexonintronprob;  // lower bound on probabilities of all exons and introns in the coding region
    int maxtracks;       // maximum reported number of transcripts at the same sequence position
    int sampleiterations;
    bool alternatives_from_sampling;
    bool alternatives_from_evidence;
    bool mea_prediction;
    bool mea_eval;
    bool needForwardTable;
    bool show_progress;
    PP::SubstateModel* profileModel;
    ContentStairs cs; // holds the local GC content class per position in the currently examined DNA
    int curGCIdx; // current index of GC content class
    list<Transcript*> *sampledTxs; // = alltranscripts stored for MultSpeciesMode
};
 
 
#endif  // _NAM_GENE_HH