website/archive/binaries/mac/src/fasta34/ag_stats.c

   1 /* this procedure implements Altschul's pre-calculated values for lambda, K */
   2
   3 /* $Name: fa_34_26_5 $ - $Id: ag_stats.c,v 1.5 2006/04/12 18:00:01 wrp Exp $ */
   4
   5 #include <stdio.h>
   6 #include <stdlib.h>
   7 #include <math.h>
   8
   9 #include "alt_parms.h"
  10
  11 static double K, Lambda, H;
  12
  13 int
  14 ag_parm(char *pam_type, int gdelval, int ggapval)
  15 {
  16   int r_v, t_gdelval, t_ggapval;
  17
  18 #ifdef OLD_FASTA_GAP
  19   t_gdelval = gdelval;
  20   t_ggapval = ggapval;
  21 #else
  22   t_gdelval = gdelval+ggapval;
  23   t_ggapval = ggapval;
  24 #endif
  25
  26   if (strcmp(pam_type,"BL50")==0 || strcmp(pam_type,"BLOSUM50")==0)
  27       r_v = look_p(bl50_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  28   else if (strcmp(pam_type,"BL62")==0 || strcmp(pam_type,"BLOSUM62")==0)
  29       r_v = look_p(bl62_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  30   else if (strcmp(pam_type,"P250")==0)
  31       r_v = look_p(p250_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  32   else if (strcmp(pam_type,"P120")==0)
  33       r_v = look_p(p120_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  34   else if (strcmp(pam_type,"MD_10")==0)
  35       r_v = look_p(md10_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  36   else if (strcmp(pam_type,"MD_20")==0)
  37       r_v = look_p(md20_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  38   else if (strcmp(pam_type,"MD_40")==0)
  39       r_v = look_p(md40_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  40   else if (strcmp(pam_type,"DNA")==0 || strcmp(pam_type,"+5/-4")==0)
  41       r_v = look_p(nt54_p,t_gdelval,t_ggapval,&K,&Lambda,&H);
  42   else r_v = 0;
  43
  44   return r_v;
  45 }
  46
  47 int
  48 look_p(struct alt_p parm[], int gap, int ext,
  49        double *K, double *Lambda, double *H)
  50 {
  51   int i;
  52
  53   gap = -gap;
  54   ext = -ext;
  55
  56   if (gap > parm[1].gap) {
  57     *K = parm[0].K;
  58     *Lambda = parm[0].Lambda;
  59     *H = parm[0].H;
  60     return 1;
  61   }
  62
  63   for (i=1; parm[i].gap > 0; i++) {
  64     if (parm[i].gap > gap) continue;
  65     else if (parm[i].gap == gap && parm[i].ext > ext ) continue;
  66     else if (parm[i].gap == gap && parm[i].ext == ext) {
  67       *K = parm[i].K;
  68       *Lambda = parm[i].Lambda;
  69       *H = parm[i].H;
  70       return 1;
  71     }
  72     else break;
  73   }
  74   return 0;
  75 }
  76
  77 int E1_to_s(double e_val, int n0, int n1) {
  78   double mp, np, a_n0, a_n0f, a_n1, a_n1f, u;
  79   int score;
  80
  81   a_n0 = (double)n0;
  82   a_n0f = log(a_n0)/H;
  83
  84   a_n1 = (double)n1;
  85   a_n1f = log(a_n1)/H;
  86
  87   mp = a_n0 - a_n0f - a_n1f;
  88   np = a_n1 - a_n0f - a_n1f;
  89
  90   if (np < 1.0) np = 1.0;
  91   if (mp < 1.0) mp = 1.0;
  92
  93   /*
  94   e_val = K * np * mp * exp ( - Lambda * score);
  95   log(e_val) = log(K np mp) - Lambda * score;
  96   (log(K np mp)-log(e_val)) / Lambda = score;
  97   */
  98   score = (int)((log( K * mp * np) - log(e_val))/Lambda +0.5);
  99   if (score < 0) score = 0;
 100   return score;
 101 }
 102
 103 double s_to_E4(int score, int n0, int  n1)
 104 {
 105   double p_val;
 106   double mp, np, a_n0, a_n0f, a_n1, a_n1f, u;
 107
 108   a_n0 = (double)n0;
 109   a_n0f = log(a_n0)/H;
 110
 111   a_n1 = (double)n1;
 112   a_n1f = log(a_n1)/H;
 113
 114   mp = a_n0 - a_n0f - a_n1f;
 115   np = a_n1 - a_n0f - a_n1f;
 116
 117   if (np < 1.0) np = 1.0;
 118   if (mp < 1.0) mp = 1.0;
 119
 120   p_val = K * np * mp * exp ( - Lambda * score);
 121
 122   if (p_val > 0.01) p_val = 1.0 - exp(-p_val);
 123
 124   return p_val * 10000.0;
 125 }
 126