1 /* this procedure implements Altschul's pre-calculated values for lambda, K */
3 /* $Name: fa_34_26_5 $ - $Id: ag_stats.c,v 1.5 2006/04/12 18:00:01 wrp Exp $ */
11 static double K, Lambda, H;
14 ag_parm(char *pam_type, int gdelval, int ggapval)
16 int r_v, t_gdelval, t_ggapval;
22 t_gdelval = gdelval+ggapval;
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);
48 look_p(struct alt_p parm[], int gap, int ext,
49 double *K, double *Lambda, double *H)
56 if (gap > parm[1].gap) {
58 *Lambda = parm[0].Lambda;
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) {
68 *Lambda = parm[i].Lambda;
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;
87 mp = a_n0 - a_n0f - a_n1f;
88 np = a_n1 - a_n0f - a_n1f;
90 if (np < 1.0) np = 1.0;
91 if (mp < 1.0) mp = 1.0;
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;
98 score = (int)((log( K * mp * np) - log(e_val))/Lambda +0.5);
99 if (score < 0) score = 0;
103 double s_to_E4(int score, int n0, int n1)
106 double mp, np, a_n0, a_n0f, a_n1, a_n1f, u;
114 mp = a_n0 - a_n0f - a_n1f;
115 np = a_n1 - a_n0f - a_n1f;
117 if (np < 1.0) np = 1.0;
118 if (mp < 1.0) mp = 1.0;
120 p_val = K * np * mp * exp ( - Lambda * score);
122 if (p_val > 0.01) p_val = 1.0 - exp(-p_val);
124 return p_val * 10000.0;