geneMSA.hh

/*
 * geneMSA.hh
 *
 * License: Artistic License, see file LICENSE.TXT or 
 *          https://opensource.org/licenses/artistic-license-1.0
 */
 
#ifndef _GENEMSA
#define _GENEMSA
 
// project includes
#include "alignment.hh"
#include "exoncand.hh"
#include "orthoexon.hh"
#include "randseqaccess.hh"
#include "speciesgraph.hh"
#include "contTimeMC.hh"
#include<unordered_map>
#include<boost/functional/hash.hpp>
 
//forward declarations
 
class OrthoGraph;
string printRFC(vector<int>);
 
class GeneMSA {
public:
    GeneMSA(RandSeqAccess *rsa, Alignment *a);
    ~GeneMSA();
 
    /*
     *  get and set functions
     */
    // start and end position (0-based, inclusive) of the gene range for the given species
    // start <= end refer to the FORWARD strand of the sequence (different from class 'alignment')
    // start = end = -1 if the species is absent from the alignment
    string getSeqID(int speciesIdx);
    Strand getStrand(int speciesIdx);
    int getStart(int speciesIdx){ return starts[speciesIdx]; }
    int getEnd(int speciesIdx){ return ends[speciesIdx]; }
    list<ExonCandidate*>* getExonCands(int speciesIdx){ return exoncands.at(speciesIdx); }
    //list<OrthoExon> getOrthoExons() { return orthoExonsList; }
    Alignment * getAlignment() {return alignment;}
    vector<int> getOffsets() {return offsets;}
 
    map<string,ExonCandidate*>* getECHash(list<ExonCandidate*> *ec); // adds the keys to the map function
 
    /*
     * assmotifqthresh, assqthresh, dssqthresh are between 0 and 1 and thresholds for the inclusion of
     * acceptor/donor splice sites based on the pattern probability
     * assqthresh=0.05 means that only acceptor ss are considered
     * that have a pattern, such that 5% of true splice site patterns have lower probability.
     * The default threshold of 0 means that all splice site patterns are considered.
     */
    void createExonCands(int s, const char *dna, map<int_fast64_t, ExonCandidate*> &ecs, map<int_fast64_t, ExonCandidate*> &addECs);
    void setExonCands(vector<map<int_fast64_t, ExonCandidate*> > &ecs);
 
    void createOrthoExons(list<OrthoExon> &orthoExonsList, map<int_fast64_t, list<pair<int,ExonCandidate*> > > &alignedECs, Evo *evo, float consThres = 0.0, int minAvLen = 0);
    void printStats(); // to stdout
    void printGeneRanges();
    void printExonCands();
    void printOrthoExons(list<OrthoExon> &orthoExonsList);
    void computeOmegas(list<OrthoExon> &orthoExonsList, vector<AnnoSequence*> const &seqRanges, PhyloTree *ctree);
    void computeOmegasEff(list<OrthoExon> &orthoExonsList, vector<AnnoSequence*> const &seqRanges, PhyloTree *ctree, ofstream *codonAli);
    void computeClamsaEff(list<OrthoExon> &orthoExonsList, vector<AnnoSequence*> const &seqRanges, PhyloTree *ctree, ofstream *codonAli);
    void computeClamsa(list<OrthoExon> &orthoExonsList, vector<AnnoSequence*> const &seqRanges, PhyloTree *ctree, ofstream *codonAli);
    vector<string> pruneToBV(vector<string> *cs, bit_vector bv); // prune codon strings to bit_vector
    vector<int> pruneToBV(vector<int> *rfc, bit_vector bv); // prune RFC to bit_vector
    double omegaForCodonTuple(vector<double> *loglik);
    void printOmegaForCodon(string outdir);
    void printCumOmega();
    void comparativeSignalScoring(list<OrthoExon> &orthoExonsList);
    // Charlotte Janas playground
    LocusTree *constructTree(); // creates, stores are returns the locus tree for the sequence tuple
 
    // calculate a columnwise conservation score and output it (for each species) in wiggle format
    void calcConsScore(list<OrthoExon> &orthoExonsList, vector<AnnoSequence*> const &seqRanges, string outdir);
    double calcColumnScore(int a, int c, int t, int g); // input: number of a,c,t and g's in one alignment column 
    void consToWig(vector<double> &consScore, string outdir);
 
    // static functions
    static void setTree(PhyloTree *t){tree = t;}
    static void setCodonEvo(CodonEvo *c){ codonevo = c; }
    static void setCodonEvoDiscr(CodonEvoDiscr *c){ codonevodiscr = c; }
    static int numSpecies(){ return tree->numSpecies(); }
    static void openOutputFiles(string outdir);
    static void closeOutputFiles();
 
    // static data members
    static int padding; // add this many bases to the region before and after the aligned region
    static int orthoExonID; // ID for exons of different species which are orthologous to each other
    static int geneRangeID; // stores an ID for the possible gene ranges of the different species which belong together
    static vector<int> exonCandID; // IDs for exon candidates of different species
    static unordered_map< bit_vector, PhyloTree*, boost::hash<bit_vector>> topologies;
    // pointers to the output files
    static vector< ofstream* > exonCands_outfiles, orthoExons_outfiles, geneRanges_outfiles_bed, geneRanges_outfiles_gff, omega_outfiles; 
    // map that stored all codon combinations on which fitch and pruning algorithm already have been run (for calculation of omega and number of substitutions)
    static map<vector<string>, pair<vector<double>, int> > computedCumValues;
 
    void printSingleOrthoExon(OrthoExon const &oe, bool files = true);
    void collect_features(int species, list<OrthoExon> *hects, SpeciesGraph *speciesgraph);
 
    void getAllOEMsas(int species, list<OrthoExon> *hects, unordered_map<string,int> *ref_class,
                         vector<AnnoSequence*> const &seqRanges);
 
    StringAlignment getMsa(OrthoExon const &oe, vector<AnnoSequence*> const &seqRanges, size_t flanking = 0);
private:
    vector<string> getCodonAlignment(OrthoExon const &oe, vector<AnnoSequence*> const &seqRanges,
                                     const vector<vector<fragment>::const_iterator > &froms,
                                     map<unsigned, vector<int> > *alignedCodons = NULL,
                                     bool generateString = true,
                                     ofstream *codonAli = NULL);
    vector<string> getCodonAlignment2(OrthoExon const &oe, vector<AnnoSequence*> const &seqRanges,
                                     const vector<vector<fragment>::const_iterator > &froms);
    void cutIncompleteCodons(OrthoExon &oe);
    cumValues* findCumValues(bit_vector bv, vector<int> rfc);
    static PhyloTree *tree;
    LocusTree *ltree;
    static CodonEvo *codonevo;
    static CodonEvoDiscr *codonevodiscr;
    vector<int> starts, ends; // gene ranges for each species
    vector<int> offsets; // this many bases are upstream from the region
    RandSeqAccess *rsa;
    Alignment* alignment;            // alignment of regions which possibly belong to a gene
    vector< list<ExonCandidate*>* > exoncands;  // exon candidates found in the different species in a gene segment
    // store cumulative omega values for every reading frame combination and every bitvector that exist 
    unordered_map<bit_vector, vector<pair<vector<int>, cumValues> >, boost::hash<bit_vector> > cumOmega; 
    map<bit_vector, map<vector<int>, vector<double> > > codonOmega; // deprecated
 
    #ifdef TESTING
    /*
    *   added by Giovanna Migliorelli 14.05.2020 
    *   code responsible for serialization
    *   aknowledgement : https://www.boost.org/doc/libs/1_70_0/libs/serialization/doc/tutorial.html
    */
 
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive & ar, const unsigned int version)
    {
        ar & starts;
        ar & alignment;
 
    }
 
    /*  added by Giovanna Migliorelli 14.05.2020
    *   default consttruct to allow for serialization
    */
 
    GeneMSA(){}; 
 
    #endif
};
 
 
 
#endif  // _GENEMSA