1 /* this procedure implements Altschul's pre-calculated values for lambda, K */
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <math.h>
6
7 #include "alt_parms.h"
8
9 static int stats_set=0;
10 static double K, Lambda, H;
11
12 int look_p(struct alt_p *parm, int gap, int ext,
13 double *K, double *Lambda, double *H);
14
15 int
set_stats(double lambda_v,double k_v)16 set_stats(double lambda_v, double k_v) {
17
18 if (lambda_v > 0.0) {Lambda=lambda_v;}
19 else return 0;
20 if (k_v > 0.0 && k_v < 1.0) {K = k_v;}
21 else return 0;
22 H = -1.0;
23 stats_set = 1;
24 return 1;
25 }
26
27 int
init_stats(char * pam_type,int * gdelval,int * ggapval,int del_set)28 init_stats(char *pam_type, int *gdelval, int *ggapval, int del_set)
29 {
30 int r_v;
31
32 struct alt_p *matrix_p;
33
34 if (stats_set) return 1;
35
36 if (strcmp(pam_type,"BL50")==0) {
37 if (!del_set) {*gdelval= -14; *ggapval=-2;}
38 matrix_p = &bl50_p[0];
39 }
40 else if (strcmp(pam_type,"BL62")==0) {
41 if (!del_set) {*gdelval= -8; *ggapval=-2;}
42 matrix_p = &bl62_p[0];
43 }
44 else if (strcmp(pam_type,"BL80")==0) {
45 if (!del_set) {*gdelval= -12; *ggapval=-2;}
46 matrix_p = &bl80_p[0];
47 }
48 else if (strcmp(pam_type,"PAM250")==0) {
49 if (!del_set) {*gdelval= -12; *ggapval=-2;}
50 matrix_p = &p250_p[0];
51 }
52 else if (strcmp(pam_type,"PAM120")==0) {
53 if (!del_set) {*gdelval= -20; *ggapval=-3;}
54 matrix_p = &p120_p[0];
55 }
56 else if (strcmp(pam_type,"MD_10")==0) {
57 if (!del_set) {*gdelval= -27; *ggapval=-4;}
58 matrix_p = &md10_p[0];
59 }
60 else if (strcmp(pam_type,"MD10")==0) {
61 if (!del_set) {*gdelval= -27; *ggapval=-4;}
62 matrix_p = &md10_p[0];
63 }
64 else if (strcmp(pam_type,"MD_20")==0) {
65 if (!del_set) {*gdelval= -26; *ggapval=-4;}
66 matrix_p = &md20_p[0];
67 }
68 else if (strcmp(pam_type,"MD20")==0) {
69 if (!del_set) {*gdelval= -26; *ggapval=-4;}
70 matrix_p = &md20_p[0];
71 }
72 else if (strcmp(pam_type,"MD_40")==0) {
73 if (!del_set) {*gdelval= -25; *ggapval=-4;}
74 matrix_p = &md40_p[0];
75 }
76 else if (strcmp(pam_type,"MD40")==0) {
77 if (!del_set) {*gdelval= -25; *ggapval=-4;}
78 matrix_p = &md40_p[0];
79 }
80 else if (strncmp(pam_type,"DNA",3)==0) {
81 if (!del_set) {*gdelval= -16; *ggapval=-4;}
82 matrix_p = &nt54_p[0];
83 }
84 else if (strncmp(pam_type,"DNA32",5)==0) {
85 if (!del_set) {*gdelval= -16; *ggapval=-4;}
86 matrix_p = &nt32_p[0];
87 }
88 else if (strncmp(pam_type,"DNA13",5)==0) {
89 if (!del_set) {*gdelval= -16; *ggapval=-4;}
90 matrix_p = &nt13_p[0];
91 }
92 else matrix_p = NULL;
93
94 if (matrix_p != NULL) {
95 #ifdef GAP_OPEN
96 if (del_set) { *gdelval += *ggapval;}
97 #endif
98 r_v = stats_set = look_p(matrix_p,*gdelval,*ggapval,&K,&Lambda,&H);
99 #ifdef GAP_OPEN
100 *gdelval -= *ggapval;
101 #endif
102 }
103 return r_v;
104 }
105
106 int
look_p(struct alt_p * parm,int gap,int ext,double * K,double * Lambda,double * H)107 look_p(struct alt_p *parm, int gap, int ext,
108 double *K, double *Lambda, double *H)
109 {
110 int i;
111
112 gap = -gap;
113 ext = -ext;
114
115 if (gap > parm[1].gap) {
116 *K = parm[0].K;
117 *Lambda = parm[0].Lambda;
118 *H = parm[0].H;
119 return 1;
120 }
121
122 for (i=1; parm[i].gap > 0; i++) {
123 if (parm[i].gap > gap) continue;
124 else if (parm[i].gap == gap && parm[i].ext > ext ) continue;
125 else if (parm[i].gap==0 || (parm[i].gap == gap && parm[i].ext == ext)) {
126 *K = parm[i].K;
127 *Lambda = parm[i].Lambda;
128 *H = parm[i].H;
129 return 1;
130 }
131 else break;
132 }
133 return 0;
134 }
135
E1_to_s(double e_val,int n0,int n1)136 int E1_to_s(double e_val, int n0, int n1) {
137 double mp, np, a_n0, a_n0f, a_n1, a_n1f, u;
138
139 if (!stats_set) return -1;
140
141 a_n0 = (double)n0;
142
143 a_n1 = (double)n1;
144
145 if (H > 0.0 ) {
146 a_n0f = log(a_n0)/H;
147 a_n1f = log(a_n1)/H;
148 }
149 else a_n0f = a_n1f = 0.0;
150
151 mp = a_n0 - a_n0f - a_n1f;
152 np = a_n1 - a_n0f - a_n1f;
153
154 if (np < 1.0) np = 1.0;
155 if (mp < 1.0) mp = 1.0;
156
157 /*
158 e_val = K * np * mp * exp ( - Lambda * score);
159 log(e_val) = log(K np mp) - Lambda * score;
160 (log(K np mp)-log(e_val)) / Lambda = score;
161 */
162 return (int)((log( K * mp * np) - log(e_val))/Lambda +0.5);
163 }
164
165 /* calculate E()-value from score, lengths, database size */
s_to_E4(int score,int n0,int n1,int z_size)166 double s_to_E4(int score, int n0, int n1, int z_size)
167 {
168 double p_val;
169 double mp, np, a_n0, a_n0f, a_n1, a_n1f, u;
170
171 if (!stats_set) return -1.0;
172
173 a_n0 = (double)n0;
174 a_n1 = (double)n1;
175
176 if (H > 0.0) {
177 a_n0f = log(a_n0)/H;
178 a_n1f = log(a_n1)/H;
179 }
180 else a_n0f = a_n1f = 0.0;
181
182 mp = a_n0 - a_n0f - a_n1f;
183 np = a_n1 - a_n0f - a_n1f;
184
185 if (np < 1.0) np = 1.0;
186 if (mp < 1.0) mp = 1.0;
187
188 p_val = K * np * mp * exp ( - Lambda * score);
189
190 if (p_val > 0.01) p_val = 1.0 - exp(-p_val);
191
192 return p_val * (double)z_size;
193 }
194
195