c_cpp/etc/KaKs_Calculator/src/GY94.h

   1 /************************************************************
   2 * Copyright (c) 2005, BGI of Chinese Academy of Sciences
   3 * All rights reserved.
   4
   5 * Filename: GY94.h
   6 * Abstract: Declaration of GY94 class.
   7
   8 * Version: 1.0
   9 * Author: Zhang Zhang  (zhanghzhang@genomics.org.cn)
  10 * Date: Oct.2, 2005
  11
  12 * Note: Source codes are taken from codeml.c in PAML.
  13
  14   References:
  15   Goldman N, Yang Z  (1994)  A codon-based model of nucleotide
  16   substitution for protein-coding DNA sequences. Mol. Biol.
  17   Evol. 11:725-736.
  18 *************************************************************/
  19
  20 #if !defined(GY94_H)
  21 #define  GY94_H
  22
  23 #include "base.h"
  24
  25 using namespace std;
  26
  27 class GY94: public Base {
  28
  29 public:
  30         GY94();
  31         GY94(string model);
  32         ~GY94();
  33
  34         /* Main function */
  35         string Run(const char *seq1, const char *seq2);
  36
  37 protected:
  38         /* Preprocess for calculating Ka&Ks */
  39         int preProcess(const char* seq1, const char* seq2);
  40         /* Parse substitution rates according to the given model */
  41         int parseSubRates(string model, double kappa[]);
  42         /* Construct two array according to genetic code */
  43         int setmark_61_64 (void);
  44         int PatternWeight();
  45         /* Encode two compared sequneces */
  46         void EncodeSeqs (void);
  47         /* Calculate Ka&Ks using ML */
  48         int PairwiseCodon (double space[]);
  49         /* Get codons' frequencies */
  50         int GetCodonFreqs(double pi[]);
  51         /* Construct transition probability matrix (64*64) */
  52         int EigenQc (int getstats, double blength, double *S, double *dS, double *dN, double Root[], double U[], double V[], double kappa[], double omega, double Q[]);
  53         /* Return maximum-likelihood score */
  54         double lfun2dSdN (double x[], int np);
  55         /* Main fuctiion for GY method */
  56         int ming2 (double *f, double x[], double xb[][2], double space[], double e, int n);
  57
  58         int eigenQREV (double Q[], double pi[], int n, double Root[], double U[], double V[], double spacesqrtpi[]);
  59         double LineSearch2 (double *f, double x0[], double p[], double step, double limit, double e, double space[], int n);
  60         double rndu(void) ;
  61         /* x[i]=x0[i] + t*p[i] */
  62         double fun_ls(double t, double x0[], double p[], double x[], int n);
  63         double distance (double x[], double y[], int n);
  64         int transform (char *z, int ls);
  65         int gradientB (int n, double x[], double f0, double g[], double space[], int xmark[]);
  66         int H_end (double x0[], double x1[], double f0, double f1, double e1, double e2, int n);
  67         void HouseholderRealSym(double a[], int n, double d[], double e[]);
  68         int EigenTridagQLImplicit(double d[], double e[], int n, double z[]);
  69         void EigenSort(double d[], double U[], int n);
  70         int eigenRealSym(double A[], int n, double Root[], double work[]);
  71         void FreeMemPUVR(void);
  72
  73 public:
  74         struct common_info {
  75                 char *z[2];     //sequences
  76                 int ns; //ns=2
  77                 int ls; //sequence's length
  78                 int ngene, npatt;
  79                 int icode;      //number of genetic code
  80                 int ncode;      //number of non-stop codon
  81                 int np, nkappa, sspace;
  82                 double fpatt[CODON*CODON];
  83                 double *space;
  84                 double kappa;   //transition/transversion
  85                 double omega;   //Ka/Ks
  86                 double pi[CODON];       //codons' frequencies
  87                 double KAPPA[8];
  88         }  com;
  89
  90 protected:
  91         double PMat[CODON*CODON],U[CODON*CODON],V[CODON*CODON],Root[CODON*CODON];
  92         int Nsensecodon, FROM61[CODON], FROM64[CODON];
  93         int Iround;
  94         unsigned int w_rndu;//=123456757;
  95         double SIZEp;
  96         double Small_Diff;
  97         string str1, str2;      //a pair of sequences
  98 };
  99
 100 #endif