c_cpp/etc/KaKs_Calculator/src/YN00.h

   1 /**********************************************************
   2 * Copyright (c) 2005, BGI of Chinese Academy of Sciences
   3 * All rights reserved.
   4
   5 * Filename: YN00.h
   6 * Abstract: Declaration of method YN00.
   7
   8 * Version: 1.0
   9 * Author: Zhang Zhang  (zhanghzhang@genomics.org.cn)
  10 * Date: Feb.2, 2005
  11
  12 * Note: Source codes are taken from yn00.c in PAML.
  13
  14   References:
  15   Yang Z, Nielsen R  (2000)  Estimating Synonymous and
  16   Nonsynonymous Substitution Rates Under Realistic
  17   Evolutionary Models. Mol Biol Evol 17:32-43.
  18 **********************************************************/
  19
  20 #if !defined(YN00_H)
  21 #define  YN00_H
  22
  23 #define CODONFREQ 12
  24 #define min(a,b) ((a)<(b)?(a):(b))
  25 #define max(a,b) ((a)>(b)?(a):(b))
  26 #define square(a) ((a)*(a))
  27 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
  28
  29 #include "base.h"
  30
  31 using namespace std;
  32
  33 class YN00: public Base {
  34
  35 public:
  36         YN00();
  37
  38         /* Main function of calculating kaks */
  39         string Run(string seq1, string seq2);
  40
  41 //protected:
  42         /* Get A,C,G,T's frequency between pair sequences: f12pos[], pi[], pi_sqrt[]  */
  43         void getFreqency(const string seq1, const string seq2);
  44         /* Get the k(transition/transversion) */
  45         virtual int GetKappa(const string seq1, const string seq2);
  46         /* Use the HKY85 Model to correct for multiple substitutions */
  47         virtual int DistanceF84(double n, double P, double Q, double pi4[],double &k_HKY, double &t, double &SEt);
  48         /* Calculate the ka,ks */
  49         virtual int DistanceYN00(const string seq1, const string seq2, double &dS,double &dN, double &SEdS, double &SEdN);
  50         /* Count synonymous and nonsynonmous sites: S, N */
  51         virtual int CountSites(const string z, double &Stot,double &Ntot,double fbS[],double fbN[]);
  52         /* Calculate the transition probability matrix using 'kappa' and 'omega' */
  53         virtual int GetPMatCodon(double P[], double kappa, double omega);
  54         /* Count synonymous and nonsynonmous differences: Sd, Nd */
  55         virtual int CountDiffs(const string seq1, const string seq2, double &Sdts,double &Sdtv,double &Ndts, double &Ndtv,double PMatrix[]);
  56
  57         //The following is for calculation of transition probability matrix by Taylor equation
  58         int eigenQREV (double Q[], double pi[], double pi_sqrt[], int n, int npi0, double Root[], double U[], double V[]);
  59         int PMatUVRoot (double P[], double t, int n, double U[], double V[], double Root[]);
  60         int eigenRealSym(double A[], int n, double Root[], double work[]);
  61         void EigenSort(double d[], double U[], int n);
  62         int EigenTridagQLImplicit(double d[], double e[], int n, double z[]);
  63         void HouseholderRealSym(double a[], int n, double d[], double e[]);
  64
  65 public:
  66         double omega;   //Ka/Ks
  67
  68 protected:
  69         /*        T     C A G
  70                 1 *     * * *
  71                 2 *     * * *
  72                 3 *     * * *
  73         f12pos: 3*4 matrix*/
  74         /* Probability of A,C,G,T at three positions */
  75         double f12pos[CODONFREQ];
  76         /* Probability of 64 codons */
  77         double pi[CODON];
  78         double pi_sqrt[CODON];
  79         /* Whether iteration calculating ka/ks */
  80         int iteration;
  81         /* The number of which pi[] is zero */
  82         int npi0;
  83
  84 };
  85
  86 #endif
  87