pp_simscore.hh

/*
 * Name:    pp_simscore.hh
 * 
 * License: Artistic License, see file LICENSE.TXT or 
 *          https://opensource.org/licenses/artistic-license-1.0
 * 
 * Project: similarity-score algorithm for block-profile and protein sequence
 * with intron information
 * Author:  Lars Gabriel
 *
 */
 
#ifndef __PP_SIMSCORE_HH
#define __PP_SIMSCORE_HH
 
#include <iostream>
#include <fstream>
#include <limits>
#include <math.h>
#include <string>
#include "properties.hh"
#include "pp_profile.hh"
 
 namespace SS {
    struct Cell {
        double score;
        std::vector < std::tuple <int, int, char > > prev;
    };
 
    class Row {
        private:
        int len;
        public:
        Cell* row;
        int position;
        Row (int length, int position);
        Cell& operator[] (int n);
        int length();
    };
 
    class SimilarityMatrix {
        private:
        std::vector<Row> matrix;
        public:
        Row& operator[] (int n);
        void addRow (int l, int p);
        int length ();
        bool empty ();
        ~SimilarityMatrix();
    };
 
    class ProteinSequence {
        private:
        char* sequence = NULL;
        std::map < int, int > introns;
        std::vector < int > num_prev_introns;
        int len;
        char* name;
        public:
        ProteinSequence (const char* fileName);
        ~ProteinSequence ();
        char& operator[] (int n);
        int length();
        int intronAt(int i);
        int intronsInRange (int start, int end);
    };
 
    struct PrflAlignment {
        //amino acid sequence that has the highest likelihood to belong
        //to the profile
        std::string argmax_aminoacids;
        //for each column the frequency of the aligned amino acid from
        //the protein sequence in the profile
        std::vector < double > match_prob;
    };
 
    struct AlignmentElement {
        int protSeq_pos;
        int block_pos;
        int col_pos;
        int length;
        char type;
    };
 
    struct BacktrackAlignDB {
        int i;
        int j;
        int j_prev;
        int block_count;
        int col_count;
        int col_posStart;
        AlignmentElement align_element;
        std::vector <AlignmentElement > alignment;
    };
 
    class SimilarityScore {
        private:
        //Similarity matrix:
        //one column represents one entry of the sequence
        //one row represents an entry of the blockprofile
        SimilarityMatrix S;
        //data structure for the protein sequence P
        ProteinSequence* seq;
        //data structure for the block profile B
        PP::Profile* prfl;
        //parameter for the gap costs
        double gap_cost_inter;
        double gap_cost_intra;
        double gap_cost_intron;
        //parameter for the intron weights
        double intron_weight_inter;
        double intron_weight_intra;
        //parameter for the pseudocount
        //the pseudocount is added to a relative frequency (v/w) of an
        //intron position with (v+epsi_intron)/(w+epsi_intron+epsi_noIntron)
        double epsi_intron;
        double epsi_noIntron;
 
        //number of protein sequences used to create the
        //intron profile of the block profile
        int intronPrfl_noProt;
 
        //background frequency of an introns at an intrablock intron position,
        //estimated from the intron positions of 15799 protein sequences
        const double intronIntraBlock_bfreq = 1.706e-3;
        //background frequency of the number of introns in
        //an intrablock intron position, estimated from the intron positions
        //of 15799 protein sequences in inter-block section
        //the length an inter-block section that was used for this estimation
        //is the mean of the inter-block distance interval
        const double intronAvgInterBlock_bfreq = 9.599e-3;
 
        //computes the intra-block intron score for a match of intron positions
        //with f residual nucleotides, if an intron of P is aligned to
        //block column s of block k
        //if f=-1: ths score of a mismatch in intron position
        //for f in {0,1,2} is computed
        double intraBlock_iscore (int k, int s, int f, int intron_frame);
        //Poisson distribution, which is used as the background distribution
        //of k intron in an inter-block section
        //lambda is estimated for I=[d_min, d_max] as
        //3 * ((d_min+d_max)/2 +1) * (intronAvgInterBlock_bfreq/intronPrfl_noProt)
        double PoiDist (int k, double lambda);
        //method for comparison, whether a score is the new best score
        //for a cell
        //and stores the information of the current best score for backtracking
        double compareScores(double new_score, double old_score, int current_i,\
            int current_j, int prev_i, int prev_j, char score_type);
        //calculates the intron score
        double intronScore (double q, double q_b);
 
        public:
        SimilarityScore (double g, double b, double g_i, double iw1,\
            double iw2, double e_i, double e_n);
        ~SimilarityScore ();
        //prints the inter block distance intervals of B to stdout
        void printInterBlock();
        bool SimilarityMatrix_empty();
        //prints an optimal alignment as the coordinates of connected
        //alignment sections of the same alignment type
        //Format: List of AlignmentElement
        void printAlignmentDB();
        //prints a list of a translations from the indice of a block to
        //the number of the block in the .prfl file.
        void printBlockPos ();
        //prints the average of the argmax of the block columns for
        //the complete profile
        pair<double, double> avgArgMaxProb();
        //read files:
        //sequence input file has to be in FASTA
        //intron positions are optional as '[Introns]' section
        //in the file format (see description of class ProteinSequence)
        //blockprofile input file has to be in .prfl
        //intron information are optional in the file format
        //(see README file)
        void readFiles (const char* seqFileName, const char* profileFileName);
 
        //algorithm to fill the similarity matrix and compute
        //the final similarity score
        void fillSimilarityMatrix ();
 
        //traces back the alignments corresponding to the final score
        //alignment are stored as list of AlignmentElements
        void backtrackingDB (int number_Alignments);
 
        //print final similarity score
        double score();
 
        //final alignments of the sequence and the profile
        std::vector <std::vector <AlignmentElement > > alignments;
 
        //datastructure for a readable alignment, which is used
        //for printing an alignment
        std::vector < pair <std::string, PrflAlignment > > alignmentsReadable;
        //prints similarity matrix and alignments with cout
        void printSimilarityMatrix ();
 
        //print readable alignment to stdout
        //FORMAT: - Alignment representation of P (AminoAcid, gap symbol
        //        or number of amino acids in inter-block)
        //        - Alignemnt representation of argmax of B
        //        (argmax AminoAcid for aligned block column,
        //        gap symbol or inter-block length)
        //        - frequency of amino acid of P in aligned block column of B,
        //        if alignment type is a match
        void printAlignment ();
    };
}
//namespace SS
#endif